The first two chunks of this r markdown file after the r setup allow for plot zooming, but it also means that the html file must be opened in a browser to view the document properly. When it knits in RStudio the preview will appear empty but the html when opened in a browser will have all the info and you can click on each plot to Zoom in on it.

Before you begin

Notes

A few notes about this script.

If you are running this with the 2022-2023 data make sure you download the whole (OSM_2022-2023 GitHub repository)[https://github.com/ACMElabUvic/OSM_2022-2023] from the ACMElabUvic GitHub. This will ensure you have all the files, data, and proper folder structure you will need to run this code and associated analyses.

Also make sure you open RStudio through the R project (OSM_2022-2023.Rproj) this will automatically set your working directory to the correct place (wherever you saved the repository) and ensure you don’t have to change the file paths for some of the data.

Lastly, if you are looking to adapt this code for a future year of data, you will want to ensure you have run the 2_ACME_landscape_covariate_exploration_script.Rmd with your data as there is much data formatting, cleaning, and restructuring that has to be done before this code will work. Helpful note: The files are numbered in the order they are used for this analysis.

If you have question please email the most recent author, currently

Marissa A. Dyck
Postdoctoral research fellow
University of Victoria
School of Environmental Studies
Email: marissadyck17@gmail.com

(update/add authors as needed)

Install packages

If you don’t already have the following packages installed, use the code below to install them.

install.packages('tidyverse') 
install.packages('ggpubr')
install.packages('corrplot')
install.packages('Hmisc')
install.packages('glmmTMB')
install.packages('MuMIn')
install.packages('TMB', type = 'source')
install.packages('rphylopic')
install.packages('broom')

Load libraries

Then load the packages to your library.

library(tidyverse) # data tidying, visualization, and much more; this will load all tidyverse packages, can see complete list using tidyverse_packages()
library(ggpubr) # make modificaions to plot for publication (arrange plots)
library(PerformanceAnalytics)    #Used to generate a correlation plot
library(Hmisc) # used to generate histograms for all variables in data frame
library(glmmTMB)      #Constructing GLMMs
## Warning in checkMatrixPackageVersion(): Package version inconsistency detected.
## TMB was built with Matrix version 1.4.1
## Current Matrix version is 1.5.3
## Please re-install 'TMB' from source using install.packages('TMB', type = 'source') or ask CRAN for a binary version of 'TMB' matching CRAN's 'Matrix' package
library(MuMIn) # for model selection
library(rphylopic) # add animal silhouettes to graphs
library(broom) # extracting odds ratios in a tidy format

Summary statistics

Load detection data

Read in saved and cleaned detection data for the 6 LUs from 2021-2022 and 2022-2023 e.g., the 1_ACME_camera_script_9-2-2024.R.

These are two separate files from the two different fiscal years, so they need to be imported and then merged into one data file for plotting. Since both are stored in the data/processed/ folder we can read them both in as a list with purrr.

# detection data
# read in saved and cleaned detection data from the ACME_camera_script_9-2-2024.R 
detections <- file.path('data/processed',
                        
                        c('OSM_ind_det_2021.csv',
                          'OSM_ind_det_2022.csv')) %>% 
  
 map(~.x %>%
        read_csv(.) %>% 
       
       # ensure the array, site, species, and event_id read in as factors
       mutate_if(is.character,
                 as.factor)) %>% 
  
  # give names to each data frame in list
  purrr::set_names('dets_2021',
                   'dets_2022') # R doesn't like when they are just numbers, you can make it work but it's annoying to call the data frame later so I've called them dets_year
## Rows: 6696 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): array, site, species, event_id
## dbl  (3): month, year, timediff
## dttm (1): datetime
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
## Rows: 14063 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (4): array, site, species, event_id
## dbl  (3): month, year, timediff
## dttm (1): datetime
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
str(detections)
## List of 2
##  $ dets_2021: tibble [6,696 × 8] (S3: tbl_df/tbl/data.frame)
##   ..$ array   : Factor w/ 2 levels "LU2","LU3": 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ site    : Factor w/ 78 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ species : Factor w/ 35 levels "Arctic hare",..: 28 28 35 35 35 35 35 35 35 35 ...
##   ..$ datetime: POSIXct[1:6696], format: "2021-07-15 11:40:07" "2022-02-06 15:11:06" ...
##   ..$ month   : num [1:6696] 7 2 7 8 8 8 8 8 8 9 ...
##   ..$ year    : num [1:6696] 2021 2022 2021 2021 2021 ...
##   ..$ timediff: num [1:6696] NA 296839 NA 28519 13230 ...
##   ..$ event_id: Factor w/ 6696 levels "E1000000","E1000001",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ dets_2022: tibble [14,063 × 8] (S3: tbl_df/tbl/data.frame)
##   ..$ array   : Factor w/ 4 levels "LU01","LU13",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ site    : Factor w/ 155 levels "LU01_06","LU01_10",..: 1 1 1 1 1 1 1 1 1 1 ...
##   ..$ species : Factor w/ 39 levels "ATVer","Beaver",..: 31 31 38 38 38 38 38 38 38 38 ...
##   ..$ datetime: POSIXct[1:14063], format: "2022-06-17 10:01:52" "2023-09-10 12:51:15" ...
##   ..$ month   : num [1:14063] 6 9 6 7 7 7 8 8 8 8 ...
##   ..$ year    : num [1:14063] 2022 2023 2022 2022 2022 ...
##   ..$ timediff: num [1:14063] NA 648166 NA 31847 21429 ...
##   ..$ event_id: Factor w/ 14063 levels "E0","E1","E10",..: 1 2 5176 6287 7398 8509 9620 10731 11842 12953 ...

Merge data

Now we need to merge these two files to make plotting easier.

They have the same columns so we could just the base R rbind() function, but in case there are differences in the columns in the future, let’s use the cleaner dplyr::bind_rows() function

# Join two years of detection data
detections_merged <- dplyr::bind_rows(detections$dets_2021,
                                      detections$dets_2022)

# check new data
str(detections_merged)
## tibble [20,759 × 8] (S3: tbl_df/tbl/data.frame)
##  $ array   : Factor w/ 6 levels "LU2","LU3","LU01",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ site    : Factor w/ 233 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ species : Factor w/ 42 levels "Arctic hare",..: 28 28 35 35 35 35 35 35 35 35 ...
##  $ datetime: POSIXct[1:20759], format: "2021-07-15 11:40:07" "2022-02-06 15:11:06" ...
##  $ month   : num [1:20759] 7 2 7 8 8 8 8 8 8 9 ...
##  $ year    : num [1:20759] 2021 2022 2021 2021 2021 ...
##  $ timediff: num [1:20759] NA 296839 NA 28519 13230 ...
##  $ event_id: Factor w/ 20759 levels "E1000000","E1000001",..: 1 2 3 4 5 6 7 8 9 10 ...

Data checks

Let’s check to make sure the data looks good after merging

Arrays and sites

Let’s check that there are the correct number of levels for array and sites, there should be 6 arrays and 233 sites (155 from 2022-2023 and 78 from 2021-2022)

str(detections_merged)
## tibble [20,759 × 8] (S3: tbl_df/tbl/data.frame)
##  $ array   : Factor w/ 6 levels "LU2","LU3","LU01",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ site    : Factor w/ 233 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ species : Factor w/ 42 levels "Arctic hare",..: 28 28 35 35 35 35 35 35 35 35 ...
##  $ datetime: POSIXct[1:20759], format: "2021-07-15 11:40:07" "2022-02-06 15:11:06" ...
##  $ month   : num [1:20759] 7 2 7 8 8 8 8 8 8 9 ...
##  $ year    : num [1:20759] 2021 2022 2021 2021 2021 ...
##  $ timediff: num [1:20759] NA 296839 NA 28519 13230 ...
##  $ event_id: Factor w/ 20759 levels "E1000000","E1000001",..: 1 2 3 4 5 6 7 8 9 10 ...
# check the levels for array and sites
levels(detections_merged$array)
## [1] "LU2"  "LU3"  "LU01" "LU13" "LU15" "LU21"
levels(detections_merged$site)
##   [1] "LU2_03"   "LU2_05"   "LU2_100"  "LU2_101"  "LU2_106"  "LU2_110" 
##   [7] "LU2_119"  "LU2_123"  "LU2_13"   "LU2_137"  "LU2_143"  "LU2_15"  
##  [13] "LU2_153"  "LU2_155"  "LU2_17"   "LU2_18"   "LU2_20"   "LU2_21"  
##  [19] "LU2_24"   "LU2_25"   "LU2_26"   "LU2_27"   "LU2_29"   "LU2_30"  
##  [25] "LU2_31"   "LU2_32"   "LU2_33"   "LU2_34"   "LU2_36"   "LU2_38"  
##  [31] "LU2_39"   "LU2_40"   "LU2_42"   "LU2_44"   "LU2_46"   "LU2_47"  
##  [37] "LU2_50"   "LU2_51"   "LU2_52"   "LU2_54"   "LU2_56"   "LU2_57"  
##  [43] "LU3_05"   "LU3_07"   "LU3_10"   "LU3_102"  "LU3_103"  "LU3_111" 
##  [49] "LU3_126"  "LU3_129"  "LU3_13"   "LU3_131"  "LU3_137"  "LU3_14"  
##  [55] "LU3_147"  "LU3_15"   "LU3_17"   "LU3_18"   "LU3_19"   "LU3_20"  
##  [61] "LU3_21"   "LU3_25"   "LU3_27"   "LU3_28"   "LU3_32"   "LU3_36"  
##  [67] "LU3_38"   "LU3_39"   "LU3_40"   "LU3_41"   "LU3_44"   "LU3_45"  
##  [73] "LU3_46"   "LU3_48"   "LU3_49"   "LU3_50"   "LU3_51"   "LU3_52"  
##  [79] "LU01_06"  "LU01_10"  "LU01_11"  "LU01_13"  "LU01_22"  "LU01_25" 
##  [85] "LU01_27"  "LU01_30"  "LU01_32"  "LU01_36"  "LU01_40"  "LU01_41" 
##  [91] "LU01_43"  "LU01_44"  "LU01_45"  "LU01_46"  "LU01_47"  "LU01_48" 
##  [97] "LU01_60"  "LU01_63"  "LU01_64"  "LU01_66"  "LU01_67"  "LU01_70" 
## [103] "LU01_71"  "LU01_72"  "LU01_73"  "LU01_74"  "LU01_75"  "LU01_76" 
## [109] "LU01_77"  "LU01_78"  "LU01_79"  "LU01_80"  "LU01_82"  "LU01_83" 
## [115] "LU01_84"  "LU01_85"  "LU01_86"  "LU13_03"  "LU13_05"  "LU13_06" 
## [121] "LU13_08"  "LU13_11"  "LU13_12"  "LU13_128" "LU13_13"  "LU13_131"
## [127] "LU13_14"  "LU13_15"  "LU13_16"  "LU13_17"  "LU13_18"  "LU13_19" 
## [133] "LU13_20"  "LU13_21"  "LU13_22"  "LU13_26"  "LU13_27"  "LU13_30" 
## [139] "LU13_32"  "LU13_33"  "LU13_34"  "LU13_35"  "LU13_36"  "LU13_37" 
## [145] "LU13_38"  "LU13_41"  "LU13_43"  "LU13_45"  "LU13_47"  "LU13_49" 
## [151] "LU13_51"  "LU13_52"  "LU13_53"  "LU13_55"  "LU13_56"  "LU13_57" 
## [157] "LU13_59"  "LU13_70"  "LU15_01"  "LU15_02"  "LU15_03"  "LU15_04" 
## [163] "LU15_07"  "LU15_08"  "LU15_09"  "LU15_10"  "LU15_11"  "LU15_12" 
## [169] "LU15_14"  "LU15_15"  "LU15_16"  "LU15_17"  "LU15_18"  "LU15_19" 
## [175] "LU15_20"  "LU15_21"  "LU15_22"  "LU15_23"  "LU15_24"  "LU15_25" 
## [181] "LU15_26"  "LU15_27"  "LU15_28"  "LU15_29"  "LU15_30"  "LU15_31" 
## [187] "LU15_32"  "LU15_34"  "LU15_36"  "LU15_37"  "LU15_40"  "LU15_41" 
## [193] "LU15_43"  "LU15_44"  "LU15_46"  "LU15_58"  "LU15_61"  "LU21_06" 
## [199] "LU21_09"  "LU21_10"  "LU21_100" "LU21_105" "LU21_106" "LU21_107"
## [205] "LU21_109" "LU21_114" "LU21_116" "LU21_119" "LU21_122" "LU21_126"
## [211] "LU21_14"  "LU21_153" "LU21_16"  "LU21_164" "LU21_21"  "LU21_23" 
## [217] "LU21_27"  "LU21_32"  "LU21_36"  "LU21_41"  "LU21_52"  "LU21_56" 
## [223] "LU21_57"  "LU21_59"  "LU21_63"  "LU21_68"  "LU21_74"  "LU21_78" 
## [229] "LU21_82"  "LU21_871" "LU21_93"  "LU21_97"  "LU21_98"

Everything looks good

NAs

Let’s ensure no NAs were introduced where there shouldn’t be during the joining process

# check for NAs introduced during data merge
summary(detections_merged)
##   array           site                    species    
##  LU2 :4711   LU01_47:  405   White-tailed deer:6141  
##  LU3 :1985   LU01_84:  396   Snowshoe hare    :4571  
##  LU01:6283   LU01_66:  383   Red squirrel     :2201  
##  LU13:1972   LU01_27:  320   Black bear       :1997  
##  LU15:2951   LU2_39 :  310   Coyote           :1285  
##  LU21:2857   LU2_50 :  292   Moose            : 693  
##              (Other):18653   (Other)          :3871  
##     datetime                          month             year     
##  Min.   :2021-07-08 10:21:40.00   Min.   : 1.000   Min.   :2021  
##  1st Qu.:2022-04-25 01:14:32.50   1st Qu.: 5.000   1st Qu.:2022  
##  Median :2022-11-03 19:10:30.00   Median : 8.000   Median :2022  
##  Mean   :2022-10-09 18:09:20.27   Mean   : 7.204   Mean   :2022  
##  3rd Qu.:2023-05-05 01:40:45.50   3rd Qu.: 9.000   3rd Qu.:2023  
##  Max.   :2023-10-02 11:08:57.00   Max.   :12.000   Max.   :2023  
##                                                                  
##     timediff          event_id    
##  Min.   :    30   E1000000:    1  
##  1st Qu.:  1427   E1000001:    1  
##  Median :  5270   E1000002:    1  
##  Mean   : 30292   E1000003:    1  
##  3rd Qu.: 18443   E1000004:    1  
##  Max.   :653883   E1000005:    1  
##  NA's   :2446     (Other) :20753

THe only NAs are in the timediff column which is what we expect since any of the first observations won’t have a timediff.

Data Formatting

In order to get plots that have the same formatting as last years’ report we have to do a bit of data formatting. First we need to make sure we are including the same relevant species (some were ignored for last years’ report or grouped together)

Last years report had the following species

  • white-tailed deer
  • snowshoe hare
  • black bear
  • coyote
  • red squirrel
  • fisher
  • unknown
  • moose
  • lynx
  • spruce grouse
  • red fox
  • striped skunk
  • ruffed grouse
  • owl
  • grey wolf
  • domestic dog
  • cougar
  • raven
  • other
  • mule deer

And they grouped all humans except for staff as ‘Humans’. Let’s look at the species we have in the combined years of data and try to format it the same way

detections_merged %>% 
  
  # group by array and species
  group_by(species) %>% 
  summarise(n = n()) %>% 
  
  # have R print everything
  print(n = nrow(.))
## # A tibble: 42 × 2
##    species                 n
##    <fct>               <int>
##  1 Arctic hare             1
##  2 ATVer                  41
##  3 Beaver                  4
##  4 Black bear           1997
##  5 Caribou               115
##  6 Cougar                 37
##  7 Coyote               1285
##  8 Domestic dog           15
##  9 Elk                     1
## 10 Fisher                257
## 11 Grey jay               66
## 12 Grey wolf             224
## 13 Human                  16
## 14 Lynx                  525
## 15 Marten                220
## 16 Moose                 693
## 17 Mule deer               2
## 18 Other                  11
## 19 Other birds           219
## 20 Owl                    15
## 21 Raven                  10
## 22 Red fox               127
## 23 Red squirrel         2201
## 24 Ruffed grouse          90
## 25 Snowmobiler            16
## 26 Snowshoe hare        4571
## 27 Spruce grouse          93
## 28 Staff                 461
## 29 Striped skunk          65
## 30 Unknown               663
## 31 Unknown canid          76
## 32 Unknown deer          340
## 33 Unknown mustelid       66
## 34 Unknown ungulate       34
## 35 White-tailed deer    6141
## 36 Canada goose            4
## 37 Hunter                  1
## 38 Long-tailed weasel     17
## 39 Otter                   7
## 40 Porcupine               5
## 41 Short-tailed weasel    19
## 42 Wolverine               8

Hmmm there is one instance of an arctic hare, check that this isn’t meant to be a snowshoe hare and fix later if needed.

Now let’s create a new data frame (tibble) to work with for the OSM figure summaries specifically

I personally would lump all the unknown together and all the birds together but for the sake of consistency with last year’s figures we will remove the same entries and keep the birds separate, let’s create a vector of entries to drop

species_drop <- c('Staff',
                  'Unknown deer',
                  'Unknown ungulate',
                  'Unknown canid',
                  'Unknown mustelid',
                  'Other birds',
                  'Arctic hare')

# now we can create the new data frame with some changes consistent w/ choices made for 2021-2022
detections_merged <- detections_merged %>% 
  
  # for summarizing, lets lump all the recreational humans into "Humans"
  mutate(species = recode_factor(species, 
                                 "Snowmobiler" = "Human",
                                 "ATVer" = "Human",
                                 'Hunter' = 'Human')) %>% 
  
  # remove species we don't want to plot
  filter(!species %in% species_drop)

# look at data
str(detections_merged)
## tibble [19,562 × 8] (S3: tbl_df/tbl/data.frame)
##  $ array   : Factor w/ 6 levels "LU2","LU3","LU01",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ site    : Factor w/ 233 levels "LU2_03","LU2_05",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ species : Factor w/ 39 levels "Human","Arctic hare",..: 33 33 33 33 33 33 33 33 33 33 ...
##  $ datetime: POSIXct[1:19562], format: "2021-07-18 08:59:30" "2021-08-07 04:21:14" ...
##  $ month   : num [1:19562] 7 8 8 8 8 8 8 9 9 9 ...
##  $ year    : num [1:19562] 2021 2021 2021 2021 2021 ...
##  $ timediff: num [1:19562] NA 28519 13230 4290 7337 ...
##  $ event_id: Factor w/ 20759 levels "E1000000","E1000001",..: 3 4 5 6 7 8 9 10 11 12 ...

Subset data

We will also want to subset the data by landscape unit (LU) and generate a new data frame for each LU to use for plotting

I’m not great at writing loops, so let’s see how this shit goes… probably bad but who knows

array_frames <- list()

for (i in unique(detections_merged$array)){
  
   #Subset data based on radius
  df <- detections_merged %>%
    filter(array == i)
  
 # list of dataframes
  array_frames <- c(array_frames, list(df))
  
 
}

# inspect one data frame
print(array_frames[[1]]) # this is for LU2
## # A tibble: 4,592 × 8
##    array site   species        datetime            month  year timediff event_id
##    <fct> <fct>  <fct>          <dttm>              <dbl> <dbl>    <dbl> <fct>   
##  1 LU2   LU2_03 White-tailed … 2021-07-18 08:59:30     7  2021      NA  E1000002
##  2 LU2   LU2_03 White-tailed … 2021-08-07 04:21:14     8  2021   28519. E1000003
##  3 LU2   LU2_03 White-tailed … 2021-08-16 08:51:21     8  2021   13230. E1000004
##  4 LU2   LU2_03 White-tailed … 2021-08-19 08:21:29     8  2021    4290. E1000005
##  5 LU2   LU2_03 White-tailed … 2021-08-24 10:39:23     8  2021    7337. E1000006
##  6 LU2   LU2_03 White-tailed … 2021-08-26 09:17:02     8  2021    2797. E1000007
##  7 LU2   LU2_03 White-tailed … 2021-08-31 19:13:33     8  2021    7796  E1000008
##  8 LU2   LU2_03 White-tailed … 2021-09-10 05:03:31     9  2021   13550. E1000009
##  9 LU2   LU2_03 White-tailed … 2021-09-16 17:10:41     9  2021    9363. E1000010
## 10 LU2   LU2_03 White-tailed … 2021-09-16 19:19:24     9  2021     127  E1000011
## # ℹ 4,582 more rows

… I think this worked

Now let’s change names of list items using purrr, couldn’t figure out how to name them in the loop, you don’t necessarily need to do this because we change the names in the next section, but I like having things named

array_frames <- array_frames %>%

  purrr::set_names('LU02',
                   'LU03',
                   'LU01',
                   'LU13',
                   'LU15',
                   'LU21') 

# inspect each data frame
head(array_frames$LU01)
## # A tibble: 6 × 8
##   array site    species        datetime            month  year timediff event_id
##   <fct> <fct>   <fct>          <dttm>              <dbl> <dbl>    <dbl> <fct>   
## 1 LU01  LU01_06 White-tailed … 2022-06-18 11:09:19     6  2022      NA  E2      
## 2 LU01  LU01_06 White-tailed … 2022-07-10 13:56:10     7  2022   31847. E3      
## 3 LU01  LU01_06 White-tailed … 2022-07-25 11:04:44     7  2022   21429. E4      
## 4 LU01  LU01_06 White-tailed … 2022-07-31 06:38:06     7  2022    8356. E5      
## 5 LU01  LU01_06 White-tailed … 2022-08-01 09:45:28     8  2022    1627. E6      
## 6 LU01  LU01_06 White-tailed … 2022-08-01 15:51:01     8  2022     364. E7

Detection plots

Detection data

Now we can apply the same data formatting for each LUs’ data frame using purrr.

We want to count the number of independent detections per species per LU to use in the detection plots

# apply the same formatting to each LU data frame using purrr map
detection_data <- array_frames %>% 
  
  purrr::map(
    ~.x %>% 
      
      # group by species
      group_by(species) %>%
      
      # calculate a column with unique accounts of each species
      mutate(count = n_distinct(event_id)) %>% 
      
      # keep just the columns we need
      select(species, count) %>% 
      
      # keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting later if you don't do it ggplot will try to count and plot each row it's annoying
      distinct()) %>% 
  
  # set names of list objects
  purrr::set_names('Detections LU02',
                   'Detections LU03',
                   'Detections LU01',
                   'Detections LU13',
                   'Detections LU15',
                   'Detections LU21')

Detection plots

Now to graph independent detections for each LU using purrr, this avoids a TON of code repetition needed to plot each one individually

We use purrr::imap() instead of purrr::map() because imap maintains the variable names in our list (e.g. Detections LU01, Detections LU13, etc.) which we can then use to title each plot.

Within purrr::imap() we just paste the code we would use for a single ggplot since all the graphical elements (except the title which we change with the file name [.y]) are the same

# create object detection plots which uses the detection_data list (w/ all 4 LUs)
detection_plots <- detection_data %>% 
  
  # use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
  purrr::imap(
    ~.x %>% 
      
      # now just copy and paste the ggplot code for the detection graphs
      ggplot(.,
             aes(x = reorder(species, count), y = count)) +
      
      # plot as bar graph using geom_col so we don't have to provide a y aesthetic
      geom_col() +
      
      # switch the x and y axis
      coord_flip() +
      
      # add the number of detections at the end of each bar
      geom_text(aes(label = count),
                color = "black",
                size = 3,
                hjust = -0.3,
                vjust = 0.2) +
      
      # label x and y axis with informative titles
      labs(x = 'Species',
           y = 'Number of Independent (30 min) Detections') +
      
      # add title to plot with LU name the .y will take the name of whatever you named each list element in the detection_data list, so make sure this name is what you want on the ggtitle
      ggtitle(.y) +
      
      # set the theme
      theme_classic() +
      theme(plot.title = element_text(hjust = 0.5)))

# view plots, this will print each in it's own window so you have to scroll back in the plot viewer pane to look at each one
detection_plots
## $`Detections LU02`

## 
## $`Detections LU03`

## 
## $`Detections LU01`

## 
## $`Detections LU13`

## 
## $`Detections LU15`

## 
## $`Detections LU21`

Save detection plots

Now we want to save these plots in case we need each individual one (we will combine the detection and naive occ plots into a single figure for each LU later and use those for the OSM report, but we may want these standalone plots later so let’s save them while they are here).

We can save all the plots from the purrr iteration above using purrr::imap. imap is used instead of map because it allows us to retain the list object names (plot names) to paste as the file name with the .y command.

IMPORTANT if you are using this code for a future github repo, DO NOT use .tiff as the file extension. This will cause issues when trying to push any changes to the github repo as the files are too large to meet githubs requirements

# save plots only use if needed
purrr::imap(
  detection_plots,
  ~ggsave(.x,
             file = paste0("figures/OSM_",
                           .y,
                           '.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
          dpi = 600,
          width = 11,
          height = 9,
          units = 'in'))

Naive occupancy

Data

We also need to alter the detection data a bit to use for naive occupancy plots.

We will use the individual LU detection data like we did before and use purrr::map() to apply the dame data formatting to all 4 data frames.

Here we want to calculate the total number of sites in each LU, the number of sites each species was detected at in each LU and then use both those numbers to calculate naive occupancy for each species in each LU

# First we need to alter the data frame a bit for these plots, let's create a data frame for each LU (I couldn't figure out how to do this without assigning individual data frames for each UGH)


# apply the same formatting to each data frame using purrr
occupancy_data <- array_frames %>% 
  
  purrr::map(
    ~.x %>% 
      
      # calculate the total number of sites for each LU
      mutate(total_sites = n_distinct(site)) %>% 
      
      # group by species to calculate the number of sites each spp occurred at
      group_by(species) %>% 
  
      # add columns to count the number of sites each spp occurred at and then the naive occupancy
  reframe(count = n_distinct(site),
          naive_occ = count/total_sites,
          ind_det = n_distinct(event_id)) %>% 
  
    # keep just the columns we need
  select(species, naive_occ, ind_det) %>% 
  
    # keep only unique (distinct) rows so we should be left with one row per species, this helps with plotting
  distinct()) %>% 
  
  purrr::set_names('Naive Occupancy LU02',
                   'Naive Occupancy LU03',
                   'Naive Occupancy LU01',
                   'Naive Occupancy LU13',
                   'Naive Occupancy LU15',
                   'Naive Occupancy LU21')

Occupancy plots

Now we can graph naive occupancy for each LU using purrr, and as with the detection plots this saves a massive amount of coding using purrr to run an iteration on the data files and produce four plots at once instead of copying and pasting code for each individually

# create object occupancy_plots which uses the occupancy_data list (w/ all 4 LUs)
occupancy_plots <- occupancy_data %>% 
  
  # use imap instead of map as it allows us to use .y to paste the list element names as the plot titles later
  purrr::imap(
    ~.x %>% 

      # now just copy and paste the ggplot code for the occupancy graphs
      ggplot(.,
             aes(x = fct_reorder(species,
                                 ind_det), # this reorders the species so they match the order of the detection plot which makes it better for viewing when the plots are arranged together in 1 figure for each LU
                 y = naive_occ)) +
      
      # plot as bars using geom_col() which uses stat = 'identity', instead of geom_bar() which will count the rows in each group and plot that instead of naive occ
      geom_col() +
      
      # flip x and y axis 
      coord_flip() +
      
      # add text to end of bars that provides naive occ value
      geom_text(aes(label = round(naive_occ, 2)), 
                size = 3, 
                hjust = -0.3, 
                vjust = 0.2) +
      
      # relabel x and y axis and title
      labs(x = 'Species',
           y = 'Proportion of Sites With At Least One Detection') +
      
      # set plot title using .y (name of list object)
      ggtitle(.y) +
      
      # set. theme elements
      theme_classic()+
      theme(plot.title = element_text(hjust = 0.5)))

# view plots
occupancy_plots
## $`Naive Occupancy LU02`

## 
## $`Naive Occupancy LU03`

## 
## $`Naive Occupancy LU01`

## 
## $`Naive Occupancy LU13`

## 
## $`Naive Occupancy LU15`

## 
## $`Naive Occupancy LU21`

Save occupancy plots

As with the detection plots, we might want these individual plots later for something so we can use purrr::imap() to save them to the figures folder

Again avoid using the .tiff extension in github

# save plots 
purrr::imap(
  occupancy_plots,
  ~ggsave(.x,
          file = paste0("figures/OSM_",
                        .y,
                        '.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
          dpi = 600,
          width = 11,
          height = 9,
          units = 'in'))

Final combined plots for report

The previous year’s report had a figure for each LU with the detections plot on the top and the occupancy plot on the bottom so we will recreate these for this year using ggarrange().

Unfortunately I could not figure out how to do this in purrr to reduce coding but luckily it isn’t too much repetition

# not sure I know how to do the following section in purrr just yet, but we've saved a ton of coding so far and it doesn't take much to arrange each of these individually

# LU2
LU02_det_occ_plots <- ggarrange(detection_plots$`Detections LU02`, 
                               occupancy_plots$`Naive Occupancy LU02`,
                               labels = c("A", "B"),
                               nrow = 2)

# view plot
LU02_det_occ_plots

# LU3
LU03_det_occ_plots <- ggarrange(detection_plots$`Detections LU03`, 
                               occupancy_plots$`Naive Occupancy LU03`,
                               labels = c("A", "B"),
                               nrow = 2)

# view plot
LU03_det_occ_plots

# LU1

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU01_det_occ_plots <- ggarrange(detection_plots$`Detections LU01`, 
                               occupancy_plots$`Naive Occupancy LU01`,
                               labels = c("A", "B"),
                               nrow = 2)

# view plot
LU01_det_occ_plots

# LU13

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU13_det_occ_plots <- ggarrange(detection_plots$`Detections LU13`, 
                                occupancy_plots$`Naive Occupancy LU13`,
                               labels = c("A", "B"),
                               nrow = 2)

# view plot
LU13_det_occ_plots

# LU15

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU15_det_occ_plots <- ggarrange(detection_plots$`Detections LU15`, 
                                occupancy_plots$`Naive Occupancy LU15`,
                                labels = c("A", "B"),
                                nrow = 2)

# view plot
LU15_det_occ_plots

# LU21

# arrange the plots so each LU has a figure with detections on top and naive occ on bottom
LU21_det_occ_plots <- ggarrange(detection_plots$`Detections LU21`, 
                                occupancy_plots$`Naive Occupancy LU21`,
                                labels = c("A", "B"),
                                nrow = 2)

# view plot
LU21_det_occ_plots

We can however, save all the figures again using purrr

# save all figures at once using purrr
final_det_occ_plots <- list(LU02_det_occ_plots,
                            LU03_det_occ_plots,
                            LU01_det_occ_plots,
                            LU13_det_occ_plots,
                            LU15_det_occ_plots,
                            LU21_det_occ_plots) %>% 
  

  purrr::set_names('LU02_det_occ_plots',
                   'LU03_det_occ_plots',
                   'LU01_det_occ_plots',
                   'LU13_det_occ_plots',
                   'LU15_det_occ_plots',
                   'LU21_det_occ_plots') %>% 
  
  purrr::imap(
    ~ggsave(.x,
            file = paste0("figures/OSM_",
                          .y,
                          '.jpg'), # avoid using .tiff extension in the github repo, those files are too large to push to origin
            dpi = 600,
            width = 12,
            height = 15,
            units = 'in'))

Finish with detection data

Before proceeding let’s clear the objects currently in our environment since we don’t need them for the analysis

rm(list = ls(all.names = TRUE))

Analysis prep

Now we can start the analysis prep.

First we need to read in the proportional detection (response metrics) and covariate (explanatory metrics) data files for all 6 LUs (fiscal years 2021-2022 and 2022-2023)

Response metrics

Read in data files

We haven’t joined the response metrics from 2021-2022 and 2022-2023 in any previous scripts yet, so let’s read those in using purrr and then join them and take a look at the data to make sure it looks good.

# response metric (proportional detections from the from the ACME_camera_script_9-2-2024.R or .Rmd)

prop_detections <-  file.path('data/processed',
                         
                         c('OSM_proportional_detections_2021.csv',
                           'OSM_proportional_detections_2022.csv')) %>% 
  
  map(~.x %>%
        read_csv(.,
                 
                  # set the column types to read in correctly
                 col_types = cols(site = col_factor(),
                                  .default = col_number()))) %>% 
  
  # give names to each data frame in list
  purrr::set_names('prop_dets_2021',
                   'prop_dets_2022') # R doesn't like when they are just numbers, you can make it work but it's annoying to call the data frame later so I've called them prop_dets_year

# check variable structure
str(prop_detections)
## List of 2
##  $ prop_dets_2021: spc_tbl_ [78 × 23] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##   ..$ site                    : Factor w/ 78 levels "LU2_03","LU2_05",..: 1 2 3 4 5 6 7 8 9 10 ...
##   ..$ black_bear              : num [1:78] 2 2 1 3 2 3 4 3 1 2 ...
##   ..$ coyote                  : num [1:78] 1 0 3 5 6 3 8 0 3 4 ...
##   ..$ fisher                  : num [1:78] 2 0 3 1 0 2 1 1 1 0 ...
##   ..$ snowshoe_hare           : num [1:78] 1 0 2 5 6 9 14 0 7 3 ...
##   ..$ white-tailed_deer       : num [1:78] 7 3 6 7 6 7 9 6 5 13 ...
##   ..$ cougar                  : num [1:78] 0 1 0 0 0 0 1 0 0 1 ...
##   ..$ lynx                    : num [1:78] 0 0 3 2 0 4 6 0 1 0 ...
##   ..$ red_fox                 : num [1:78] 0 0 3 0 0 0 0 1 0 0 ...
##   ..$ moose                   : num [1:78] 0 0 0 0 1 0 3 3 0 1 ...
##   ..$ grey_wolf               : num [1:78] 0 0 0 0 0 0 0 0 1 0 ...
##   ..$ caribou                 : num [1:78] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ absent_black_bear       : num [1:78] 3 3 4 2 3 8 7 2 4 9 ...
##   ..$ absent_coyote           : num [1:78] 6 7 4 2 1 11 6 7 4 10 ...
##   ..$ absent_fisher           : num [1:78] 5 7 4 6 7 12 13 6 6 14 ...
##   ..$ absent_snowshoe_hare    : num [1:78] 6 7 5 2 1 5 0 7 0 11 ...
##   ..$ absent_white-tailed_deer: num [1:78] 0 4 1 0 1 7 5 1 2 1 ...
##   ..$ absent_cougar           : num [1:78] 7 6 7 7 7 14 13 7 7 13 ...
##   ..$ absent_lynx             : num [1:78] 7 7 4 5 7 10 8 7 6 14 ...
##   ..$ absent_red_fox          : num [1:78] 7 7 4 7 7 14 14 6 7 14 ...
##   ..$ absent_moose            : num [1:78] 7 7 7 7 6 14 11 4 7 13 ...
##   ..$ absent_grey_wolf        : num [1:78] 7 7 7 7 7 14 14 7 6 14 ...
##   ..$ absent_caribou          : num [1:78] 7 7 7 7 7 14 14 7 7 14 ...
##   ..- attr(*, "spec")=
##   .. .. cols(
##   .. ..   .default = col_number(),
##   .. ..   site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   .. ..   black_bear = col_number(),
##   .. ..   coyote = col_number(),
##   .. ..   fisher = col_number(),
##   .. ..   snowshoe_hare = col_number(),
##   .. ..   `white-tailed_deer` = col_number(),
##   .. ..   cougar = col_number(),
##   .. ..   lynx = col_number(),
##   .. ..   red_fox = col_number(),
##   .. ..   moose = col_number(),
##   .. ..   grey_wolf = col_number(),
##   .. ..   caribou = col_number(),
##   .. ..   absent_black_bear = col_number(),
##   .. ..   absent_coyote = col_number(),
##   .. ..   absent_fisher = col_number(),
##   .. ..   absent_snowshoe_hare = col_number(),
##   .. ..   `absent_white-tailed_deer` = col_number(),
##   .. ..   absent_cougar = col_number(),
##   .. ..   absent_lynx = col_number(),
##   .. ..   absent_red_fox = col_number(),
##   .. ..   absent_moose = col_number(),
##   .. ..   absent_grey_wolf = col_number(),
##   .. ..   absent_caribou = col_number()
##   .. .. )
##   ..- attr(*, "problems")=<externalptr> 
##  $ prop_dets_2022: spc_tbl_ [154 × 25] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##   ..$ site                    : Factor w/ 154 levels "LU01_06","LU01_10",..: 1 2 3 4 5 6 7 8 9 10 ...
##   ..$ black_bear              : num [1:154] 7 3 4 7 8 9 4 5 7 7 ...
##   ..$ coyote                  : num [1:154] 4 4 8 10 11 9 11 0 9 4 ...
##   ..$ fisher                  : num [1:154] 4 3 3 3 2 1 1 2 0 3 ...
##   ..$ moose                   : num [1:154] 3 2 5 9 1 0 2 4 1 0 ...
##   ..$ snowshoe_hare           : num [1:154] 4 1 3 0 8 2 2 0 12 4 ...
##   ..$ white-tailed_deer       : num [1:154] 12 5 12 12 13 14 15 9 12 10 ...
##   ..$ cougar                  : num [1:154] 0 0 1 0 1 0 0 0 0 0 ...
##   ..$ grey_wolf               : num [1:154] 0 0 2 0 0 0 1 0 0 0 ...
##   ..$ lynx                    : num [1:154] 0 0 1 0 1 1 0 0 0 2 ...
##   ..$ red_fox                 : num [1:154] 0 0 2 0 0 0 0 0 4 0 ...
##   ..$ wolverine               : num [1:154] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ caribou                 : num [1:154] 0 0 0 0 0 0 0 0 0 0 ...
##   ..$ absent_black_bear       : num [1:154] 5 3 8 5 4 3 8 7 5 5 ...
##   ..$ absent_coyote           : num [1:154] 10 1 6 5 3 5 4 15 6 11 ...
##   ..$ absent_fisher           : num [1:154] 10 2 11 12 12 13 14 13 15 12 ...
##   ..$ absent_moose            : num [1:154] 11 3 9 6 13 14 13 11 14 15 ...
##   ..$ absent_snowshoe_hare    : num [1:154] 10 4 11 15 6 12 13 15 3 11 ...
##   ..$ absent_white-tailed_deer: num [1:154] 2 0 2 3 1 0 0 6 3 5 ...
##   ..$ absent_cougar           : num [1:154] 14 5 13 15 13 14 15 15 15 15 ...
##   ..$ absent_grey_wolf        : num [1:154] 14 5 12 15 14 14 14 15 15 15 ...
##   ..$ absent_lynx             : num [1:154] 14 5 13 15 13 13 15 15 15 13 ...
##   ..$ absent_red_fox          : num [1:154] 14 5 12 15 14 14 15 15 11 15 ...
##   ..$ absent_wolverine        : num [1:154] 14 5 14 15 14 14 15 15 15 15 ...
##   ..$ absent_caribou          : num [1:154] 14 5 14 15 14 14 15 15 15 15 ...
##   ..- attr(*, "spec")=
##   .. .. cols(
##   .. ..   .default = col_number(),
##   .. ..   site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   .. ..   black_bear = col_number(),
##   .. ..   coyote = col_number(),
##   .. ..   fisher = col_number(),
##   .. ..   moose = col_number(),
##   .. ..   snowshoe_hare = col_number(),
##   .. ..   `white-tailed_deer` = col_number(),
##   .. ..   cougar = col_number(),
##   .. ..   grey_wolf = col_number(),
##   .. ..   lynx = col_number(),
##   .. ..   red_fox = col_number(),
##   .. ..   wolverine = col_number(),
##   .. ..   caribou = col_number(),
##   .. ..   absent_black_bear = col_number(),
##   .. ..   absent_coyote = col_number(),
##   .. ..   absent_fisher = col_number(),
##   .. ..   absent_moose = col_number(),
##   .. ..   absent_snowshoe_hare = col_number(),
##   .. ..   `absent_white-tailed_deer` = col_number(),
##   .. ..   absent_cougar = col_number(),
##   .. ..   absent_grey_wolf = col_number(),
##   .. ..   absent_lynx = col_number(),
##   .. ..   absent_red_fox = col_number(),
##   .. ..   absent_wolverine = col_number(),
##   .. ..   absent_caribou = col_number()
##   .. .. )
##   ..- attr(*, "problems")=<externalptr>

Merge data files

Now we need to merge the two data files for analysis. We can do this with dplyr.

# merge the proportional detections files so there are rows for both fiscal years
prop_dets_all <- dplyr::bind_rows(prop_detections$prop_dets_2021,
                                  prop_detections$prop_dets_2022)

print(prop_dets_all)
## # A tibble: 232 × 25
##    site  black_bear coyote fisher snowshoe_hare `white-tailed_deer` cougar  lynx
##    <fct>      <dbl>  <dbl>  <dbl>         <dbl>               <dbl>  <dbl> <dbl>
##  1 LU2_…          2      1      2             1                   7      0     0
##  2 LU2_…          2      0      0             0                   3      1     0
##  3 LU2_…          1      3      3             2                   6      0     3
##  4 LU2_…          3      5      1             5                   7      0     2
##  5 LU2_…          2      6      0             6                   6      0     0
##  6 LU2_…          3      3      2             9                   7      0     4
##  7 LU2_…          4      8      1            14                   9      1     6
##  8 LU2_…          3      0      1             0                   6      0     0
##  9 LU2_…          1      3      1             7                   5      0     1
## 10 LU2_…          2      4      0             3                  13      1     0
## # ℹ 222 more rows
## # ℹ 17 more variables: red_fox <dbl>, moose <dbl>, grey_wolf <dbl>,
## #   caribou <dbl>, absent_black_bear <dbl>, absent_coyote <dbl>,
## #   absent_fisher <dbl>, absent_snowshoe_hare <dbl>,
## #   `absent_white-tailed_deer` <dbl>, absent_cougar <dbl>, absent_lynx <dbl>,
## #   absent_red_fox <dbl>, absent_moose <dbl>, absent_grey_wolf <dbl>,
## #   absent_caribou <dbl>, wolverine <dbl>, absent_wolverine <dbl>

Format data

This looks good except since there were no wolverines in the first fiscal year of monitoring (LU02 and LU03) those columns have NAs for both arrays, we want to replace those NAs with Zeros and move the wolverine column to the correct location

Let’s do that now.

prop_dets_all <- prop_dets_all %>% 
  
   # replace NAs introduced from joining data to zeros
  replace(is.na(.),
          0) %>% 
  
  # relocate wolverine column
  relocate(.,
           wolverine,
           .after = caribou)

# check data
head(prop_dets_all)
## # A tibble: 6 × 25
##   site   black_bear coyote fisher snowshoe_hare `white-tailed_deer` cougar  lynx
##   <fct>       <dbl>  <dbl>  <dbl>         <dbl>               <dbl>  <dbl> <dbl>
## 1 LU2_03          2      1      2             1                   7      0     0
## 2 LU2_05          2      0      0             0                   3      1     0
## 3 LU2_1…          1      3      3             2                   6      0     3
## 4 LU2_1…          3      5      1             5                   7      0     2
## 5 LU2_1…          2      6      0             6                   6      0     0
## 6 LU2_1…          3      3      2             9                   7      0     4
## # ℹ 17 more variables: red_fox <dbl>, moose <dbl>, grey_wolf <dbl>,
## #   caribou <dbl>, wolverine <dbl>, absent_black_bear <dbl>,
## #   absent_coyote <dbl>, absent_fisher <dbl>, absent_snowshoe_hare <dbl>,
## #   `absent_white-tailed_deer` <dbl>, absent_cougar <dbl>, absent_lynx <dbl>,
## #   absent_red_fox <dbl>, absent_moose <dbl>, absent_grey_wolf <dbl>,
## #   absent_caribou <dbl>, absent_wolverine <dbl>

Looks good!

Save data

Let’s save the merged and formatted detection data from 2021 and 2022 for future use

write_csv(prop_dets_all,
          'data/processed/OSM_proportional_detections_merged_2021_2022.csv')

Covariates

In the previous script, 2_ACME_landscape_covariate_exploration_script.Rmd we joined the two fiscal years of data and grouped the covariates for analysis and saved this data as a c.sv file, so we can read in this file now and we shouldn’t have to do any further formatting at the moment

Read data

We will check the data structure after reading in the file just to make sure everything looks good.

covariates_all <- read_csv('data/processed/OSM_covariates_grouped_2021_2022.csv',
                           
                           # set the column types to read in correctly
                           col_types = cols(array = col_factor(),
                                            camera = col_factor(),
                                            site = col_factor(),
                                            buff_dist = col_factor(),
                                            .default = col_number()))
## Warning: The following named parsers don't match the column names: camera
str(covariates_all)
## spc_tbl_ [4,660 × 19] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ array             : Factor w/ 6 levels "LU13","LU15",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ site              : Factor w/ 233 levels "LU13_18","LU13_15",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ buff_dist         : Factor w/ 20 levels "250","500","750",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ harvest           : num [1:4660] 0 0 0.687 0.337 0 ...
##  $ pipeline          : num [1:4660] 0 0.068 0 0 0.0301 ...
##  $ roads             : num [1:4660] 0 0.0174 0 0 0 ...
##  $ seismic_lines     : num [1:4660] 0 0.03277 0 0.00889 0.01144 ...
##  $ seismic_lines_3D  : num [1:4660] 0 0 0 0 0.0523 ...
##  $ trails            : num [1:4660] 0.00588 0.0028 0 0.01591 0 ...
##  $ transmission_lines: num [1:4660] 0.0642 0 0 0 0.091 ...
##  $ veg_edges         : num [1:4660] 0 0.0858 0 0 0 ...
##  $ wells             : num [1:4660] 0 0 0 0 0.0322 ...
##  $ lc_grassland      : num [1:4660] 0.193 0.348 0 0 0.178 ...
##  $ lc_coniferous     : num [1:4660] 0.456 0.358 0.186 1 0.822 ...
##  $ lc_broadleaf      : num [1:4660] 0 0 0 0 0 ...
##  $ lc_mixed          : num [1:4660] 0 0.101 0.255 0 0 ...
##  $ lc_developed      : num [1:4660] 0 0.0916 0 0 0 ...
##  $ lc_shrub          : num [1:4660] 0.316 0 0.559 0 0 ...
##  $ osm_industrial    : num [1:4660] 0.383 0.157 0 0 0 ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   .default = col_number(),
##   ..   array = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   site = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   buff_dist = col_factor(levels = NULL, ordered = FALSE, include_na = FALSE),
##   ..   harvest = col_number(),
##   ..   pipeline = col_number(),
##   ..   roads = col_number(),
##   ..   seismic_lines = col_number(),
##   ..   seismic_lines_3D = col_number(),
##   ..   trails = col_number(),
##   ..   transmission_lines = col_number(),
##   ..   veg_edges = col_number(),
##   ..   wells = col_number(),
##   ..   lc_grassland = col_number(),
##   ..   lc_coniferous = col_number(),
##   ..   lc_broadleaf = col_number(),
##   ..   lc_mixed = col_number(),
##   ..   lc_developed = col_number(),
##   ..   lc_shrub = col_number(),
##   ..   osm_industrial = col_number()
##   .. )
##  - attr(*, "problems")=<externalptr>

Everything looks good!

###Subset data by buffer

We need to subset the data so we have separate data frames for each buffer width to work with in the analysis AND to explore correlation between variables at each buffer width, as these may very with spatial scales

Let’s use a for loop to subset the data

buffer_frames <- list()

for (i in unique(covariates_all$buff_dist)){
  
  print(i)
  
  # Subset data based on radius
  df <- covariates_all %>%
    filter(buff_dist == i)
  
  # list of dataframes
  buffer_frames <-c (buffer_frames, list(df))
}
## [1] "250"
## [1] "500"
## [1] "750"
## [1] "1000"
## [1] "1250"
## [1] "1500"
## [1] "1750"
## [1] "2000"
## [1] "2250"
## [1] "2500"
## [1] "2750"
## [1] "3000"
## [1] "3250"
## [1] "3500"
## [1] "3750"
## [1] "4000"
## [1] "4250"
## [1] "4500"
## [1] "4750"
## [1] "5000"
# name list objects so we can extract names for plotting 

buffer_frames <- buffer_frames %>% 
  
  # absurdly long way to do this but for sake of time fuck it
  purrr::set_names('250 meter buffer',
                   '500 meter buffer',
                   '750 meter buffer',
                   '1000 meter buffer',
                   '1250 meter buffer',
                   '1500 meter buffer',
                   '1750 meter buffer',
                   '2000 meter buffer',
                   '2250 meter buffer',
                   '2500 meter buffer',
                   '2750 meter buffer',
                   '3000 meter buffer',
                   '3250 meter buffer',
                   '3500 meter buffer',
                   '3750 meter buffer',
                   '4000 meter buffer',
                   '4250 meter buffer',
                   '4500 meter buffer',
                   '4750 meter buffer',
                   '5000 meter buffer')

Now we have a list with data frames for each buffer width which we can work with later.

Add response metric

Now that we have the covariate data formatted we need to add the response metric (monthly proportional presence/absence) to the data frames

osm_final_df_2021_2022 <- buffer_frames %>% 
  
  purrr::map(
    ~.x %>% 
      
      left_join(prop_dets_all,
                by = 'site'))

Finish with data formatting

Let’s remove the objects we no longer need from the environment to keep our workspace clean

rm(covariates_all,
   prop_detections,
   df,
   i)

Analysis

Now we are going to run a global model which includes all HFI and LC variables that at first glance (will do a more thorough check later) seem to have enough data to include as covariates for each buffer width, and then we will compare these models see which buffer width best fit the data for each species. After that we will optimize models so they don’t includes any variables that are highly correlated.

We don’t need to do ALL the species since many don’t have enough data.

Refer to the 1_ACME_camera_script_9-2-2024.html or .Rmd the plot for proportional monthly detections should provide info on which species we have enough data for, can be found under Response metrics/3.Proportion monthly detections

A brief look at this fig indicates that we have enough for all the mammals in the prop_detections data frame except

  • cougar
  • wolverine
  • caribou??? (may have enough, may not)

Black bear

there is probably a way to shorten the following code to select particular species, I saw Andrew’s for loop in the draft script he wrote but couldn’t quite figure out how to adapt it to my purposes with the data formatted the way I have it, so I did this instead, maybe we can merge approaches later to clean this up if deemed necessary? But it certainly functions for now and is understandable… I think.

Global models

Let’s start with bears and use purrr to create a global model for every buffer distance

Recall purrr::map() is magical for iterations and will apply all the functions within the map() function to each item of the list supplied before the the map() function.

# create models for black bears at each buffer size
black_bear_mods <- osm_final_df_2021_2022 %>%

  # use purrr map to run the same functions for all buffer sizes ((all objects in list))
  purrr::map(
    ~.x %>%

      # glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
     glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))

Model selection

We will use the model.sel() function from the MuMIn package to compare the global models for each buffer width and see which buffer fits the black bear data best

model.sel(black_bear_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 250 meter buffer     -0.5998          +      0.012140           0.1221
## 1750 meter buffer    -0.5987          +      0.008720           0.3954
## 2000 meter buffer    -0.5976          +      0.005699           0.4518
## 1250 meter buffer    -0.5996          +      0.008946           0.3455
## 1500 meter buffer    -0.6002          +      0.016130           0.3382
## 1000 meter buffer    -0.5983          +     -0.001997           0.4015
## 5000 meter buffer    -0.6074          +      0.047290           0.1153
## 2250 meter buffer    -0.5960          +     -0.001043           0.4491
## 2750 meter buffer    -0.5959          +      0.005900           0.4160
## 2500 meter buffer    -0.5954          +      0.001508           0.4295
## 3750 meter buffer    -0.6021          +      0.024810           0.3098
## 4000 meter buffer    -0.6037          +      0.035760           0.3021
## 3500 meter buffer    -0.6004          +      0.023550           0.3449
## 3250 meter buffer    -0.5989          +      0.024520           0.3668
## 4750 meter buffer    -0.6059          +      0.037070           0.1895
## 4250 meter buffer    -0.6040          +      0.040600           0.2960
## 3000 meter buffer    -0.5967          +      0.011800           0.3862
## 500 meter buffer     -0.6006          +      0.029270           0.2784
## 4500 meter buffer    -0.6048          +      0.037590           0.2450
## 750 meter buffer     -0.5978          +      0.022290           0.3135
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 250 meter buffer           0.01706        -0.086020         0.067700
## 1750 meter buffer          0.28860        -0.149500        -0.001904
## 2000 meter buffer          0.37870        -0.137500         0.004301
## 1250 meter buffer          0.18040        -0.146700        -0.007598
## 1500 meter buffer          0.18560        -0.166800        -0.010420
## 1000 meter buffer          0.25690        -0.112100         0.021270
## 5000 meter buffer         -0.05003        -0.101400         0.093600
## 2250 meter buffer          0.39070        -0.153800        -0.004265
## 2750 meter buffer          0.35480        -0.169300         0.014090
## 2500 meter buffer          0.37300        -0.152200        -0.006222
## 3750 meter buffer          0.17330        -0.159600         0.034770
## 4000 meter buffer          0.15840        -0.140500         0.061200
## 3500 meter buffer          0.22620        -0.155000         0.028350
## 3250 meter buffer          0.27570        -0.156100         0.024090
## 4750 meter buffer          0.03390        -0.103500         0.088520
## 4250 meter buffer          0.14460        -0.117100         0.074410
## 3000 meter buffer          0.31280        -0.166900         0.015930
## 500 meter buffer           0.15610        -0.001471         0.108800
## 4500 meter buffer          0.08703        -0.113800         0.082690
## 750 meter buffer           0.18220        -0.048030         0.044580
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 250 meter buffer          0.054090          0.11220          0.046220
## 1750 meter buffer         0.150400          0.17960         -0.023060
## 2000 meter buffer         0.195400          0.22580         -0.030690
## 1250 meter buffer         0.133300          0.15110          0.026870
## 1500 meter buffer         0.114200          0.14190          0.004238
## 1000 meter buffer         0.184100          0.17820         -0.002263
## 5000 meter buffer         0.008365          0.04685         -0.077250
## 2250 meter buffer         0.205400          0.23040         -0.045190
## 2750 meter buffer         0.189300          0.21420         -0.071020
## 2500 meter buffer         0.197600          0.22710         -0.056540
## 3750 meter buffer         0.093200          0.14440         -0.106300
## 4000 meter buffer         0.078050          0.13870         -0.085970
## 3500 meter buffer         0.112300          0.16360         -0.106600
## 3250 meter buffer         0.133400          0.17710         -0.106000
## 4750 meter buffer         0.044460          0.08694         -0.084660
## 4250 meter buffer         0.071410          0.14580         -0.070480
## 3000 meter buffer         0.167100          0.19230         -0.090640
## 500 meter buffer          0.114000          0.14040          0.023410
## 4500 meter buffer         0.056320          0.11480         -0.080600
## 750 meter buffer          0.138500          0.11560          0.035930
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 250 meter buffer      4.783e-02     -0.104000             -0.11560
## 1750 meter buffer    -7.498e-02      0.158300             -0.11170
## 2000 meter buffer    -9.514e-02      0.176000             -0.10060
## 1250 meter buffer     9.031e-03     -0.006562             -0.13600
## 1500 meter buffer    -2.918e-02      0.073920             -0.12820
## 1000 meter buffer    -9.271e-03      0.017020             -0.12980
## 5000 meter buffer    -4.227e-01      0.005939             -0.08420
## 2250 meter buffer    -9.128e-02      0.202100             -0.09823
## 2750 meter buffer    -1.476e-01      0.273700             -0.08796
## 2500 meter buffer    -1.030e-01      0.245000             -0.09330
## 3750 meter buffer    -1.610e-01      0.117700             -0.08006
## 4000 meter buffer    -2.101e-01      0.062950             -0.06989
## 3500 meter buffer    -1.497e-01      0.174000             -0.07614
## 3250 meter buffer    -1.592e-01      0.241900             -0.07557
## 4750 meter buffer    -3.426e-01      0.037880             -0.07104
## 4250 meter buffer    -2.148e-01      0.031390             -0.06951
## 3000 meter buffer    -1.587e-01      0.273800             -0.08137
## 500 meter buffer      7.892e-05     -0.133900             -0.10010
## 4500 meter buffer    -2.709e-01      0.042870             -0.06702
## 750 meter buffer      2.574e-02     -0.074510             -0.12570
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 250 meter buffer           -0.03550       0.10260         -0.129200
## 1750 meter buffer          -0.18240       0.07213         -0.092160
## 2000 meter buffer          -0.17810       0.07224         -0.078630
## 1250 meter buffer          -0.13410       0.08228         -0.128800
## 1500 meter buffer          -0.16190       0.07507         -0.114800
## 1000 meter buffer          -0.12410       0.09118         -0.114300
## 5000 meter buffer          -0.13330       0.01758          0.224000
## 2250 meter buffer          -0.17130       0.07448         -0.093590
## 2750 meter buffer          -0.16380       0.05180         -0.067080
## 2500 meter buffer          -0.16810       0.06390         -0.086760
## 3750 meter buffer          -0.15790       0.07173         -0.004969
## 4000 meter buffer          -0.15410       0.06589          0.029670
## 3500 meter buffer          -0.15900       0.06431         -0.033120
## 3250 meter buffer          -0.15790       0.05717         -0.044450
## 4750 meter buffer          -0.14610       0.03566          0.154500
## 4250 meter buffer          -0.16110       0.05966          0.052200
## 3000 meter buffer          -0.15550       0.04647         -0.055050
## 500 meter buffer           -0.04670       0.05432         -0.104000
## 4500 meter buffer          -0.15530       0.05068          0.096340
## 750 meter buffer           -0.06456       0.04651         -0.107200
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 250 meter buffer          -0.009170       -0.0136 18 -443.782 926.8  0.00
## 1750 meter buffer         -0.047050        0.1401 18 -444.092 927.4  0.62
## 2000 meter buffer         -0.062020        0.1507 18 -444.158 927.5  0.75
## 1250 meter buffer          0.033960        0.1370 18 -444.226 927.7  0.89
## 1500 meter buffer         -0.008584        0.1543 18 -444.304 927.8  1.04
## 1000 meter buffer          0.027230        0.1196 18 -444.336 927.9  1.11
## 5000 meter buffer         -0.102300        0.4047 18 -444.816 928.8  2.07
## 2250 meter buffer         -0.058860        0.1602 18 -445.030 929.3  2.50
## 2750 meter buffer         -0.099740        0.1762 18 -445.409 930.0  3.25
## 2500 meter buffer         -0.093760        0.1587 18 -445.417 930.0  3.27
## 3750 meter buffer         -0.052440        0.2391 18 -445.548 930.3  3.53
## 4000 meter buffer         -0.053510        0.2714 18 -445.553 930.3  3.54
## 3500 meter buffer         -0.069670        0.2051 18 -445.859 930.9  4.15
## 3250 meter buffer         -0.096130        0.1884 18 -445.949 931.1  4.33
## 4750 meter buffer         -0.086190        0.3327 18 -446.118 931.4  4.67
## 4250 meter buffer         -0.062220        0.2681 18 -446.193 931.6  4.82
## 3000 meter buffer         -0.100500        0.1798 18 -446.234 931.7  4.90
## 500 meter buffer           0.024600        0.0979 18 -446.296 931.8  5.03
## 4500 meter buffer         -0.073450        0.2919 18 -446.364 931.9  5.16
## 750 meter buffer          -0.003772        0.1360 18 -447.921 935.1  8.28
##                   weight
## 250 meter buffer   0.159
## 1750 meter buffer  0.117
## 2000 meter buffer  0.109
## 1250 meter buffer  0.102
## 1500 meter buffer  0.094
## 1000 meter buffer  0.091
## 5000 meter buffer  0.057
## 2250 meter buffer  0.046
## 2750 meter buffer  0.031
## 2500 meter buffer  0.031
## 3750 meter buffer  0.027
## 4000 meter buffer  0.027
## 3500 meter buffer  0.020
## 3250 meter buffer  0.018
## 4750 meter buffer  0.015
## 4250 meter buffer  0.014
## 3000 meter buffer  0.014
## 500 meter buffer   0.013
## 4500 meter buffer  0.012
## 750 meter buffer   0.003
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

Looks like the smallest buffer (250) fits the data best for black bears.

Let’s look at this model closer

Model summary

summary(black_bear_mods$`250 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(black_bear, absent_black_bear) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    923.6    985.6   -443.8    887.6      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.05453  0.2335  
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -0.59976    0.10675  -5.618 1.93e-08 ***
## scale(harvest)             0.01214    0.05072   0.239   0.8109    
## scale(pipeline)            0.04783    0.07103   0.673   0.5007    
## scale(roads)              -0.10401    0.09070  -1.147   0.2515    
## scale(seismic_lines)      -0.03550    0.05305  -0.669   0.5034    
## scale(seismic_lines_3D)   -0.11563    0.05881  -1.966   0.0493 *  
## scale(trails)              0.10263    0.04771   2.151   0.0315 *  
## scale(transmission_lines) -0.12923    0.06773  -1.908   0.0564 .  
## scale(veg_edges)          -0.00917    0.07788  -0.118   0.9063    
## scale(wells)              -0.01360    0.05402  -0.252   0.8012    
## scale(osm_industrial)      0.04622    0.05292   0.873   0.3825    
## scale(lc_grassland)        0.06770    0.09942   0.681   0.4959    
## scale(lc_coniferous)       0.01706    0.26332   0.065   0.9483    
## scale(lc_broadleaf)        0.12207    0.18033   0.677   0.4985    
## scale(lc_mixed)            0.05409    0.15802   0.342   0.7321    
## scale(lc_developed)       -0.08602    0.09897  -0.869   0.3847    
## scale(lc_shrub)            0.11218    0.14259   0.787   0.4314    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Nothing looks fishy in the model summary for now, we will look at this more closely once we have a true top model.

Code to remove a model

At one point the 250 meter buffer was giving us issues so we had to remove it. After re-extracting the data and re-doing some data formatting this is no longer an issue but I’ve saved the code here in case it’s needed in the future

# create models for black bears at each buffer size
black_bear_mods_no250 <- osm_final_df_2021_2022 %>%

# remove 250 meter buffer width
  purrr::discard_at('250 meter buffer') %>% 

  # use purrr map to fun the same functions for all buffer sizes ((all objects in list))
  purrr::map(
    ~.x %>%

      # glmmTMB function let's us run the proportional binomial model using cbind to combine the present and absent columns for each species
     glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))

Subset models

Autocorrelation 250m

Before we can develop model subsets we need to see what variables can be included in the same model at this buffer width.

Let’s use the chart.Correlation() function in the Performance Analytics package to look at this.

buffer_frames$`250 meter buffer` %>% 
  
  select_if(is.numeric) %>% 
  
   # use chart.correlation 
      chart.Correlation(.,
                        histogram = TRUE, 
                        method = "pearson")

mtext('250 meter buffer', side = 3, line = 3)

> You can click on this fig to zoom in!

List of correlated variables:

  • pipeline & transmission_lines 0.53
  • roads & veg_edges 0.71
  • roads & lc_developed 0.57

Model list (Andrew can you review this)

Jake/Andrew can you suggest some additional subset models if needed

Let’s create another global model without these correlated variables. I’m going to select transmission_lines over pipeline because the summary from earlier showed transmission lines had larger effect on black bear presence, and I’m going to choose to keep roads instead of veg edges and the developed landcover class because we are interested in the effect of roads more than these other two variables.

# global model w/ non-correlated variables
bear_global_250 <-  glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = osm_final_df_2021_2022$`250 meter buffer`,
                   family = 'binomial')

# null model to compare
bear_null_250 <- glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~ 1 +
                        
                        # Random effect of array
                        (1|array),
                   data = osm_final_df_2021_2022$`250 meter buffer`,
                   family = 'binomial')

# second model is based on linear features providing easier movement through boreal forest
bear_linear_250 <- glmmTMB::glmmTMB(cbind(black_bear, absent_black_bear) ~
                        
                        # HFI
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        
                        # Random effect of array
                        (1|array),
                   data = osm_final_df_2021_2022$`250 meter buffer`,
                   family = 'binomial')

Model selection

Now let’s look at a model selection table with our subset models.

model.sel(bear_global_250,
          bear_null_250,
          bear_linear_250)
## Model selection table 
##                 cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## bear_linear_250    -0.5986          +                               
## bear_global_250    -0.5995          +       0.01223           0.2326
## bear_null_250      -0.5884          +                               
##                 cnd(scl(lc_cnf)) cnd(scl(lc_grs)) cnd(scl(lc_mxd))
## bear_linear_250                                                   
## bear_global_250           0.1823           0.1248           0.1455
## bear_null_250                                                     
##                 cnd(scl(lc_shr)) cnd(scl(osm_ind)) cnd(scl(rds))
## bear_linear_250                                          -0.1586
## bear_global_250           0.1891           0.04922       -0.1215
## bear_null_250                                                   
##                 cnd(scl(ssm_lns_3D)) cnd(scl(ssm_lns)) cnd(scl(trl))
## bear_linear_250              -0.1311          -0.04436       0.10250
## bear_global_250              -0.1135          -0.04271       0.09811
## bear_null_250                                                       
##                 cnd(scl(trn_lns)) cnd(scl(wll)) df   logLik  AICc delta weight
## bear_linear_250          -0.09901                7 -449.674 913.8  0.00  0.972
## bear_global_250          -0.10520      -0.01348 15 -444.383 921.0  7.14  0.027
## bear_null_250                                    2 -463.090 930.2 16.38  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

Looks like the linear feature model is best, SO FAR

Top model

250 meter buffer - linear features

Summary

summary(bear_linear_250)
##  Family: binomial  ( logit )
## Formula:          
## cbind(black_bear, absent_black_bear) ~ scale(roads) + scale(seismic_lines) +  
##     scale(seismic_lines_3D) + scale(trails) + scale(transmission_lines) +  
##     (1 | array)
## Data: osm_final_df_2021_2022$`250 meter buffer`
## 
##      AIC      BIC   logLik deviance df.resid 
##    913.3    937.5   -449.7    899.3      225 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.05426  0.2329  
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -0.59856    0.10644  -5.624 1.87e-08 ***
## scale(roads)              -0.15864    0.05694  -2.786  0.00534 ** 
## scale(seismic_lines)      -0.04436    0.05006  -0.886  0.37556    
## scale(seismic_lines_3D)   -0.13114    0.05635  -2.327  0.01995 *  
## scale(trails)              0.10247    0.04641   2.208  0.02724 *  
## scale(transmission_lines) -0.09901    0.05653  -1.751  0.07987 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Odds ratios

Let’s extract the odds ratios for the top model so we can plot them for data vis later.

bear_model_odds <- bear_linear_250 %>% 

  # extract the coefficients and upper and lower CI
      confint() %>% 
 
      # format resulting object as a tibble data frame
      as_tibble() %>%
  
      
      # add a column where we can put the feature names
      rowid_to_column() %>% 
      
      # rename the columns for plotting
      rename('lower' = `2.5 %`,
             'upper' = `97.5 %`,
             'estimate' = Estimate,
             'feature' = rowid) %>% 
      
      # rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
      mutate(feature = as.factor(feature),
             feature = recode(feature,
                              '1' = 'intercept',
                              '2' = 'roads',
                              '3' = 'seismic_lines',
                              '4' = 'seismic_lines_3D',
                              '5' = 'trails',
                              '6' = 'transmission_lines',
                              '7' = 'intercept_array')) %>% 
  
  # remove intercepts
  filter(!grepl('intercept',
                feature))

Plot odds ratios

First let’s get a silhouette for this graphy from phylopic

black_bear_img <- get_phylopic(get_uuid(name = 'Ursus americanus'))

Now let’s use ggplot to plot the odds ratios for each feature in the top model

# provide data and mapping aesthetics
ggplot(bear_model_odds, aes(x = feature, 
                          y = estimate)) +
  
   geom_errorbar(aes(ymin = lower, 
                     ymax = upper),
                width = 0.4,
                linewidth = 0.5,
                position = position_dodge(width = 0.9)) +
  
   geom_hline(yintercept = 0, linetype = "dashed") +
  
  labs(y = 'Odds ratio') +
  
  scale_x_discrete(labels = c('Roads',
                              'Seismic Lines',
                              '3D Seismic Lines',
                              'Trails',
                              'Transmission Lines')) +
  
  add_phylopic(black_bear_img, 
               x = 5.3,
               y = 0.2,
               ysize = 0.05) +
  
  theme_classic() +
  
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 90,
                                   hjust = 1),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 14))

Let’s repeat this process for each species that we have enough data for.

Caribou

We may or may not have enough data for caribou but let’s try it at least for this preliminary report

We can use the same code from black bears (above) to run global models for each buffer width except remember we want to remove 250 meters

And in the same chunk to save time let’s also run the model.sel() function

Global models

caribou_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(caribou, absent_caribou) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))
## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')

## Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
## problem; singular convergence (7). See vignette('troubleshooting'),
## help('diagnose')
# model selection
model.sel(caribou_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 3000 meter buffer     -8.327          +       0.25940           2.3660
## 3250 meter buffer     -8.413          +       0.40820           4.7700
## 2500 meter buffer     -7.922          +       0.02584          -0.2111
## 3500 meter buffer     -8.443          +       0.58490           5.7480
## 2750 meter buffer     -7.879          +       0.11420           1.0120
## 1750 meter buffer     -7.391          +      -0.05851           1.9440
## 2250 meter buffer     -7.581          +       0.06022           0.0175
## 2000 meter buffer     -7.068          +      -0.09839           0.3358
## 1500 meter buffer     -7.190          +      -0.14750           0.8697
## 3750 meter buffer     -7.991          +       0.68070           6.2590
## 1250 meter buffer    -71.530          +      -0.24830          -0.8018
## 5000 meter buffer     -7.286          +      -0.21240           6.0100
## 4000 meter buffer     -7.325          +       0.84700           6.9050
## 4750 meter buffer     -7.182          +       0.34290           7.0820
## 4500 meter buffer     -7.264          +       0.96420           7.9520
## 4250 meter buffer     -7.112          +       0.93660           7.1020
## 1000 meter buffer   -121.300          +      -0.88000          -1.4550
## 750 meter buffer     -43.960          +      -1.42900          -1.6690
## 500 meter buffer    -840.800          +      -2.41900          -1.4380
## 250 meter buffer    -206.200          +    -421.40000          -2.7170
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 3000 meter buffer          10.3800           3.4700           1.2090
## 3250 meter buffer          12.5800           3.6880           1.6250
## 2500 meter buffer           6.9770           2.5380           0.4845
## 3500 meter buffer          12.5000           3.8100           1.6340
## 2750 meter buffer           8.3090           2.9310           0.7919
## 1750 meter buffer           9.4270           2.1820           1.0480
## 2250 meter buffer           6.6580           2.1220           0.4970
## 2000 meter buffer           6.4720           1.9540           0.5902
## 1500 meter buffer           7.0280           1.3290           0.8263
## 3750 meter buffer          12.0400           3.3390           1.5220
## 1250 meter buffer           3.1020           0.1831           0.3333
## 5000 meter buffer          12.0500           2.1740           1.7240
## 4000 meter buffer          12.0400           2.2370           1.4110
## 4750 meter buffer          12.4500           1.5170           1.6060
## 4500 meter buffer          13.2500           1.4420           1.6180
## 4250 meter buffer          12.1800           1.5420           1.4130
## 1000 meter buffer           0.5239          -0.1699          -0.1059
## 750 meter buffer           -0.7779          -0.2413          -0.3634
## 500 meter buffer           -0.5981          -0.5873          -0.3747
## 250 meter buffer           -0.9887          -0.9365          -0.5740
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 3000 meter buffer          4.58300          5.10400           0.04095
## 3250 meter buffer          4.75200          6.13000          -0.07332
## 2500 meter buffer          3.41100          3.49100          -0.30210
## 3500 meter buffer          3.81300          6.00900          -0.07073
## 2750 meter buffer          3.85900          4.12200          -0.02226
## 1750 meter buffer          4.86600          4.94700          -0.49040
## 2250 meter buffer          3.07400          3.33900          -0.69810
## 2000 meter buffer          3.25200          3.35000          -0.60850
## 1500 meter buffer          3.56400          3.61300          -0.53180
## 3750 meter buffer          3.16900          5.75500          -0.02580
## 1250 meter buffer          1.25600          1.42400          -0.97300
## 5000 meter buffer          3.74900          6.21200           1.64700
## 4000 meter buffer          3.21900          5.96400          -0.05433
## 4750 meter buffer          3.45400          6.38300           1.05400
## 4500 meter buffer          3.57200          6.74800           0.16230
## 4250 meter buffer          3.30100          6.14000          -0.04194
## 1000 meter buffer         -0.05488         -0.02029          -1.44600
## 750 meter buffer          -0.75030         -0.73140          -0.97290
## 500 meter buffer          -0.52240         -0.60790          -0.24610
## 250 meter buffer          -0.67780         -0.72120          -0.08231
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 3000 meter buffer      -0.87190       1.11100             0.058070
## 3250 meter buffer      -1.06100       1.73800             0.408500
## 2500 meter buffer      -0.02702      -0.22180            -0.172100
## 3500 meter buffer      -1.31300       1.85600             0.596900
## 2750 meter buffer      -0.45420       0.30130            -0.147300
## 1750 meter buffer       0.91290      -0.64880            -0.005105
## 2250 meter buffer       0.48070      -0.10070            -0.096940
## 2000 meter buffer       0.50400      -0.32360            -0.039980
## 1500 meter buffer       0.97430      -0.10050             0.011860
## 3750 meter buffer      -1.35300       1.67200             0.751900
## 1250 meter buffer       0.59150       0.61250             0.047650
## 5000 meter buffer      -0.39520      -1.02000             0.284300
## 4000 meter buffer      -0.73410       0.96520             0.631600
## 4750 meter buffer      -0.38340      -0.41850             0.277700
## 4500 meter buffer      -0.39080       0.72670             0.413400
## 4250 meter buffer      -0.52920       0.99310             0.467600
## 1000 meter buffer       0.51580       0.77400             0.024020
## 750 meter buffer        0.42510       0.10420            -0.097430
## 500 meter buffer        0.16950       0.32310            -0.199400
## 250 meter buffer        0.07895       0.05935            -0.291100
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 3000 meter buffer           0.36310     -0.202000           -1.7600
## 3250 meter buffer           0.49490     -0.291400           -2.0620
## 2500 meter buffer           0.43350      0.006096           -1.0940
## 3500 meter buffer           0.68140     -0.388600           -2.2220
## 2750 meter buffer           0.38360     -0.083580           -1.2860
## 1750 meter buffer          -0.14810      0.010880           -1.7330
## 2250 meter buffer           0.45640      0.147300           -1.1670
## 2000 meter buffer          -0.03146      0.064910           -1.4170
## 1500 meter buffer          -0.08270      0.060250           -2.3370
## 3750 meter buffer           0.68810     -0.366600           -1.8540
## 1250 meter buffer           0.22600      0.245600         -146.9000
## 5000 meter buffer           0.33590     -0.334100           -0.8979
## 4000 meter buffer           0.52040     -0.324100           -1.0870
## 4750 meter buffer           0.54200     -0.384300           -0.5912
## 4500 meter buffer           0.57490     -0.332600           -0.5100
## 4250 meter buffer           0.51920     -0.347700           -0.6046
## 1000 meter buffer           0.20020      0.245300         -274.0000
## 750 meter buffer            0.25540      0.159600          -94.4300
## 500 meter buffer            0.27440      0.041970        -2395.0000
## 250 meter buffer            0.35630     -0.269500         -202.0000
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 3000 meter buffer          -2.92600       1.53200 18 -106.914 253.0  0.00
## 3250 meter buffer          -3.27800       1.24400 18 -107.147 253.5  0.47
## 2500 meter buffer          -1.75900       1.61500 18 -107.456 254.1  1.09
## 3500 meter buffer          -3.68400       1.25900 18 -107.532 254.3  1.24
## 2750 meter buffer          -2.21100       1.55700 18 -108.019 255.2  2.21
## 1750 meter buffer          -0.85030       1.28600 18 -108.491 256.2  3.15
## 2250 meter buffer          -1.39300       1.16900 18 -109.426 258.1  5.02
## 2000 meter buffer          -1.30400       1.29900 18 -110.065 259.3  6.30
## 1500 meter buffer          -0.88710       0.95770 18 -110.314 259.8  6.80
## 3750 meter buffer          -2.87300       1.12000 18 -110.797 260.8  7.77
## 1250 meter buffer          -0.16140       0.38550 18 -112.914 265.0 12.00
## 5000 meter buffer           0.05748       0.25730 18 -114.018 267.2 14.21
## 4000 meter buffer          -0.44760       0.61300 18 -114.115 267.4 14.40
## 4750 meter buffer           0.38050       0.20930 18 -115.073 269.4 16.32
## 4500 meter buffer           0.54040       0.09117 18 -115.505 270.2 17.18
## 4250 meter buffer           0.26280       0.27810 18 -115.784 270.8 17.74
## 1000 meter buffer          -0.10720       0.30250 18 -116.045 271.3 18.26
## 750 meter buffer           -0.16640       0.45710 18 -121.221 281.7 28.61
## 500 meter buffer            0.17530       0.38700 18 -127.857 294.9 41.89
## 250 meter buffer            0.02991       0.27320 18 -129.695 298.6 45.56
##                   weight
## 3000 meter buffer  0.275
## 3250 meter buffer  0.218
## 2500 meter buffer  0.160
## 3500 meter buffer  0.148
## 2750 meter buffer  0.091
## 1750 meter buffer  0.057
## 2250 meter buffer  0.022
## 2000 meter buffer  0.012
## 1500 meter buffer  0.009
## 3750 meter buffer  0.006
## 1250 meter buffer  0.001
## 5000 meter buffer  0.000
## 4000 meter buffer  0.000
## 4750 meter buffer  0.000
## 4500 meter buffer  0.000
## 4250 meter buffer  0.000
## 1000 meter buffer  0.000
## 750 meter buffer   0.000
## 500 meter buffer   0.000
## 250 meter buffer   0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

We get a warning that there are some model convergence problems, I expect this is because we don’t have enough data for caribou but I don’t have time to dig into this now so we will investigate more closely for final analysis.

For now let’s examine the top buffer width and see if the model seems reasonable.

For caribou 3000m buffer is top model for now

Model summary 3000m

Let’s take a closer look at the top model summary

summary(caribou_mods$`3000 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(caribou, absent_caribou) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    249.8    311.9   -106.9    213.8      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.823    0.9072  
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -8.32675    1.10886  -7.509 5.95e-14 ***
## scale(harvest)             0.25940    0.59626   0.435  0.66352    
## scale(pipeline)           -0.87191    0.46951  -1.857  0.06330 .  
## scale(roads)               1.11082    1.11248   0.999  0.31803    
## scale(seismic_lines)       0.36310    0.27271   1.331  0.18305    
## scale(seismic_lines_3D)    0.05807    0.24641   0.236  0.81370    
## scale(trails)             -0.20199    0.27759  -0.728  0.46683    
## scale(transmission_lines) -1.76031    0.96766  -1.819  0.06889 .  
## scale(veg_edges)          -2.92633    0.97157  -3.012  0.00260 ** 
## scale(wells)               1.53157    0.51256   2.988  0.00281 ** 
## scale(osm_industrial)      0.04095    0.66097   0.062  0.95060    
## scale(lc_grassland)        1.20931    0.71225   1.698  0.08953 .  
## scale(lc_coniferous)      10.37851    4.95897   2.093  0.03636 *  
## scale(lc_broadleaf)        2.36607    3.34699   0.707  0.47961    
## scale(lc_mixed)            4.58332    2.47197   1.854  0.06372 .  
## scale(lc_developed)        3.47022    1.10097   3.152  0.00162 ** 
## scale(lc_shrub)            5.10412    2.65456   1.923  0.05451 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There’s nothing that catches my eye immediately as being sus about this particular model so it may not have been one with convergence issues. We will keep it in report for now

Subset Models

Add later if deemed necessary.

Top Model

Coyote

Global models

coyote_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(coyote, absent_coyote) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(coyote_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 2250 meter buffer     -1.395          +      0.042500          -0.6109
## 2500 meter buffer     -1.395          +      0.042930          -0.6368
## 2000 meter buffer     -1.392          +      0.041820          -0.5804
## 2750 meter buffer     -1.392          +      0.032190          -0.6176
## 3000 meter buffer     -1.392          +      0.027240          -0.6311
## 1750 meter buffer     -1.392          +      0.056400          -0.5381
## 3250 meter buffer     -1.391          +      0.016470          -0.6451
## 1500 meter buffer     -1.392          +      0.064670          -0.5189
## 4750 meter buffer     -1.392          +     -0.059700          -0.6479
## 3500 meter buffer     -1.389          +      0.003136          -0.6677
## 5000 meter buffer     -1.392          +     -0.065830          -0.6166
## 4500 meter buffer     -1.390          +     -0.051360          -0.6518
## 1250 meter buffer     -1.396          +      0.046620          -0.5118
## 3750 meter buffer     -1.388          +      0.007226          -0.6813
## 4250 meter buffer     -1.388          +     -0.028700          -0.6442
## 4000 meter buffer     -1.387          +     -0.007103          -0.6526
## 1000 meter buffer     -1.401          +      0.039700          -0.4491
## 750 meter buffer      -1.407          +      0.071470          -0.4037
## 500 meter buffer      -1.405          +      0.026170          -0.2139
## 250 meter buffer      -1.401          +     -0.070110          -0.2204
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 2250 meter buffer          -0.9458          0.15070         -0.23690
## 2500 meter buffer          -0.9737          0.11920         -0.23530
## 2000 meter buffer          -0.9085          0.17180         -0.23420
## 2750 meter buffer          -0.9469          0.09757         -0.20270
## 3000 meter buffer          -0.9893          0.09057         -0.19160
## 1750 meter buffer          -0.8653          0.19850         -0.20550
## 3250 meter buffer          -0.9989          0.09923         -0.19700
## 1500 meter buffer          -0.8416          0.20920         -0.16630
## 4750 meter buffer          -0.9970          0.13270         -0.30750
## 3500 meter buffer          -1.0740          0.11500         -0.23170
## 5000 meter buffer          -0.9541          0.17180         -0.30810
## 4500 meter buffer          -1.0420          0.14040         -0.31040
## 1250 meter buffer          -0.7872          0.23960         -0.15930
## 3750 meter buffer          -1.1050          0.12710         -0.24990
## 4250 meter buffer          -1.0320          0.12570         -0.28310
## 4000 meter buffer          -1.0800          0.10770         -0.25700
## 1000 meter buffer          -0.6724          0.29370         -0.13670
## 750 meter buffer           -0.5802          0.31140         -0.10340
## 500 meter buffer           -0.3256          0.26610         -0.01793
## 250 meter buffer           -0.3305          0.26450         -0.04606
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 2250 meter buffer         -0.52840         -0.48660           0.15630
## 2500 meter buffer         -0.52510         -0.50680           0.12900
## 2000 meter buffer         -0.51090         -0.45900           0.19800
## 2750 meter buffer         -0.49540         -0.48320           0.12130
## 3000 meter buffer         -0.50530         -0.48900           0.11380
## 1750 meter buffer         -0.50070         -0.41700           0.21220
## 3250 meter buffer         -0.49670         -0.48300           0.11090
## 1500 meter buffer         -0.50420         -0.39330           0.24070
## 4750 meter buffer         -0.40570         -0.43400           0.23000
## 3500 meter buffer         -0.51840         -0.50790           0.08270
## 5000 meter buffer         -0.38280         -0.41390           0.25370
## 4500 meter buffer         -0.43660         -0.44890           0.20870
## 1250 meter buffer         -0.48650         -0.37750           0.23710
## 3750 meter buffer         -0.53730         -0.50340           0.09756
## 4250 meter buffer         -0.45700         -0.44640           0.19250
## 4000 meter buffer         -0.51430         -0.47740           0.13990
## 1000 meter buffer         -0.41680         -0.28340           0.26170
## 750 meter buffer          -0.34620         -0.23400           0.29250
## 500 meter buffer          -0.16310         -0.07484           0.26540
## 250 meter buffer          -0.08185         -0.05793           0.20090
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 2250 meter buffer      -0.19420       0.79880             -0.05081
## 2500 meter buffer      -0.20700       0.77430             -0.04427
## 2000 meter buffer      -0.12810       0.72600             -0.06077
## 2750 meter buffer      -0.23170       0.74130             -0.03743
## 3000 meter buffer      -0.20270       0.63160             -0.02955
## 1750 meter buffer      -0.08028       0.61160             -0.08205
## 3250 meter buffer      -0.15920       0.50920             -0.02971
## 1500 meter buffer      -0.04714       0.45540             -0.10910
## 4750 meter buffer       0.02514       0.11170             -0.09480
## 3500 meter buffer      -0.11030       0.42650             -0.03691
## 5000 meter buffer       0.03415       0.07007             -0.09585
## 4500 meter buffer       0.05460       0.12340             -0.09086
## 1250 meter buffer      -0.00669       0.39200             -0.14390
## 3750 meter buffer      -0.08029       0.33360             -0.04971
## 4250 meter buffer       0.02218       0.15910             -0.08018
## 4000 meter buffer      -0.02849       0.23380             -0.06312
## 1000 meter buffer       0.05279       0.26340             -0.14670
## 750 meter buffer        0.09082       0.13220             -0.14820
## 500 meter buffer        0.02605       0.01657             -0.11120
## 250 meter buffer        0.11220      -0.14140             -0.03538
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 2250 meter buffer            0.3687     0.0278600           0.15570
## 2500 meter buffer            0.3613     0.0129400           0.14920
## 2000 meter buffer            0.3399     0.0342100           0.08978
## 2750 meter buffer            0.3269     0.0012820           0.12700
## 3000 meter buffer            0.3186     0.0017140           0.09048
## 1750 meter buffer            0.3195     0.0472500           0.07428
## 3250 meter buffer            0.3124     0.0038310           0.08569
## 1500 meter buffer            0.3121     0.0418300           0.05657
## 4750 meter buffer            0.3328    -0.0236500           0.08922
## 3500 meter buffer            0.3291     0.0088240           0.07006
## 5000 meter buffer            0.3392    -0.0148300           0.09121
## 4500 meter buffer            0.3408    -0.0149500           0.05294
## 1250 meter buffer            0.3072     0.0224100          -0.01659
## 3750 meter buffer            0.3300     0.0007599           0.07923
## 4250 meter buffer            0.3238    -0.0146400           0.06037
## 4000 meter buffer            0.3275    -0.0057330           0.06099
## 1000 meter buffer            0.2722     0.0020390          -0.04874
## 750 meter buffer             0.2598    -0.0421500           0.02405
## 500 meter buffer             0.1991    -0.0910900           0.08040
## 250 meter buffer             0.1448    -0.0115700           0.02799
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik   AICc delta
## 2250 meter buffer          -0.43050    -0.0005123 18 -446.617  932.4  0.00
## 2500 meter buffer          -0.38210     0.0412500 18 -447.591  934.4  1.95
## 2000 meter buffer          -0.43650     0.0037450 18 -449.597  938.4  5.96
## 2750 meter buffer          -0.30720     0.0468400 18 -450.058  939.3  6.88
## 3000 meter buffer          -0.20420     0.0654700 18 -451.851  942.9 10.47
## 1750 meter buffer          -0.38830     0.0145800 18 -452.179  943.6 11.12
## 3250 meter buffer          -0.11390     0.0657700 18 -454.425  948.1 15.62
## 1500 meter buffer          -0.29670     0.0411900 18 -455.179  949.6 17.12
## 4750 meter buffer           0.12040     0.1108000 18 -455.299  949.8 17.36
## 3500 meter buffer          -0.04928     0.0933100 18 -455.355  949.9 17.47
## 5000 meter buffer           0.10430     0.1126000 18 -455.580  950.4 17.93
## 4500 meter buffer           0.10020     0.1188000 18 -456.109  951.4 18.98
## 1250 meter buffer          -0.25020     0.0282600 18 -456.244  951.7 19.25
## 3750 meter buffer          -0.01394     0.1228000 18 -457.485  954.2 21.74
## 4250 meter buffer           0.08758     0.1119000 18 -458.373  956.0 23.51
## 4000 meter buffer           0.06694     0.1143000 18 -459.143  957.5 25.05
## 1000 meter buffer          -0.24580     0.0890100 18 -460.179  959.6 27.12
## 750 meter buffer           -0.26040     0.1760000 18 -463.775  966.8 34.31
## 500 meter buffer           -0.24080     0.2046000 18 -478.613  996.4 63.99
## 250 meter buffer           -0.16590     0.0354700 18 -493.595 1026.4 93.96
##                   weight
## 2250 meter buffer  0.680
## 2500 meter buffer  0.257
## 2000 meter buffer  0.035
## 2750 meter buffer  0.022
## 3000 meter buffer  0.004
## 1750 meter buffer  0.003
## 3250 meter buffer  0.000
## 1500 meter buffer  0.000
## 4750 meter buffer  0.000
## 3500 meter buffer  0.000
## 5000 meter buffer  0.000
## 4500 meter buffer  0.000
## 1250 meter buffer  0.000
## 3750 meter buffer  0.000
## 4250 meter buffer  0.000
## 4000 meter buffer  0.000
## 1000 meter buffer  0.000
## 750 meter buffer   0.000
## 500 meter buffer   0.000
## 250 meter buffer   0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

for coyote top model appears to be 2250 m

Model summary 2250m

Let’s get the model summary for this model

summary(coyote_mods$`2250 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(coyote, absent_coyote) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    929.2    991.3   -446.6    893.2      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.04572  0.2138  
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.3953463  0.1033031 -13.507  < 2e-16 ***
## scale(harvest)             0.0424989  0.0787950   0.539 0.589638    
## scale(pipeline)           -0.1942326  0.1112808  -1.745 0.080911 .  
## scale(roads)               0.7988319  0.1681427   4.751 2.02e-06 ***
## scale(seismic_lines)       0.3687452  0.0685662   5.378 7.53e-08 ***
## scale(seismic_lines_3D)   -0.0508061  0.0633194  -0.802 0.422335    
## scale(trails)              0.0278614  0.0639774   0.435 0.663208    
## scale(transmission_lines)  0.1557429  0.0947854   1.643 0.100360    
## scale(veg_edges)          -0.4305253  0.1445842  -2.978 0.002904 ** 
## scale(wells)              -0.0005123  0.1048918  -0.005 0.996103    
## scale(osm_industrial)      0.1562622  0.0724134   2.158 0.030934 *  
## scale(lc_grassland)       -0.2368575  0.0814280  -2.909 0.003628 ** 
## scale(lc_coniferous)      -0.9457625  0.2442679  -3.872 0.000108 ***
## scale(lc_broadleaf)       -0.6109090  0.1683292  -3.629 0.000284 ***
## scale(lc_mixed)           -0.5283815  0.1511077  -3.497 0.000471 ***
## scale(lc_developed)        0.1506740  0.1159080   1.300 0.193620    
## scale(lc_shrub)           -0.4865559  0.1392808  -3.493 0.000477 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Everything looks good here

Subset models

Add later if deemed necessary.

Top model

Fisher

Global models

fisher_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(fisher, absent_fisher) ~
                        
                         # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))
# model selection
model.sel(fisher_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 1000 meter buffer     -3.103          +       -0.2130           0.7220
## 2250 meter buffer     -3.032          +       -0.1094           0.4305
## 2500 meter buffer     -3.015          +       -0.1106           0.4519
## 750 meter buffer      -3.073          +       -0.2519           0.9443
## 500 meter buffer      -3.061          +       -0.3083           0.8742
## 2000 meter buffer     -3.030          +       -0.1302           0.4442
## 2750 meter buffer     -3.009          +       -0.1313           0.5002
## 3000 meter buffer     -3.010          +       -0.1380           0.5026
## 1250 meter buffer     -3.057          +       -0.2260           0.5567
## 3250 meter buffer     -3.014          +       -0.1605           0.4713
## 4250 meter buffer     -3.049          +       -0.1686           0.2283
## 4500 meter buffer     -3.043          +       -0.1821           0.2328
## 1500 meter buffer     -3.021          +       -0.2114           0.4784
## 4750 meter buffer     -3.054          +       -0.2156           0.2438
## 1750 meter buffer     -3.017          +       -0.1744           0.4529
## 4000 meter buffer     -3.044          +       -0.1668           0.2790
## 250 meter buffer      -3.050          +       -0.3459           0.0651
## 3500 meter buffer     -3.016          +       -0.1812           0.4238
## 3750 meter buffer     -3.034          +       -0.1816           0.3237
## 5000 meter buffer     -3.035          +       -0.2367           0.2430
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 1000 meter buffer           0.9855           0.4920          0.17270
## 2250 meter buffer           0.6056           0.3171         -0.06029
## 2500 meter buffer           0.6796           0.3220         -0.13240
## 750 meter buffer            1.2770           0.5516          0.25020
## 500 meter buffer            1.1280           0.5176          0.28000
## 2000 meter buffer           0.6137           0.3267         -0.05118
## 2750 meter buffer           0.8144           0.4164         -0.16220
## 3000 meter buffer           0.8642           0.5128         -0.20340
## 1250 meter buffer           0.7391           0.3690          0.10130
## 3250 meter buffer           0.8923           0.5791         -0.21840
## 4250 meter buffer           0.6988           0.8591         -0.14580
## 4500 meter buffer           0.6921           0.9013         -0.11230
## 1500 meter buffer           0.6358           0.2761          0.03435
## 4750 meter buffer           0.7411           0.8550         -0.06847
## 1750 meter buffer           0.6202           0.2558         -0.02444
## 4000 meter buffer           0.7327           0.7886         -0.12320
## 250 meter buffer           -0.0501           0.3446          0.03051
## 3500 meter buffer           0.8674           0.6212         -0.18970
## 3750 meter buffer           0.7937           0.7089         -0.14090
## 5000 meter buffer           0.6989           0.8561         -0.05248
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 1000 meter buffer          0.43740       -0.1534000          -0.08007
## 2250 meter buffer          0.24660       -0.2819000          -0.17220
## 2500 meter buffer          0.30300       -0.1841000          -0.17230
## 750 meter buffer           0.65500        0.0708200          -0.09660
## 500 meter buffer           0.58170        0.0360900          -0.21500
## 2000 meter buffer          0.25140       -0.2750000          -0.15340
## 2750 meter buffer          0.38860       -0.0505900          -0.11080
## 3000 meter buffer          0.40150        0.0002221          -0.08809
## 1250 meter buffer          0.32890       -0.2144000          -0.13390
## 3250 meter buffer          0.43260        0.0062080          -0.15360
## 4250 meter buffer          0.37530       -0.1945000          -0.38390
## 4500 meter buffer          0.40890       -0.1884000          -0.37390
## 1500 meter buffer          0.27800       -0.1981000          -0.14170
## 4750 meter buffer          0.46060       -0.1798000          -0.47500
## 1750 meter buffer          0.24740       -0.2467000          -0.13110
## 4000 meter buffer          0.35640       -0.1960000          -0.30560
## 250 meter buffer          -0.00188       -0.4457000          -0.25220
## 3500 meter buffer          0.42070       -0.0375800          -0.22720
## 3750 meter buffer          0.38160       -0.1261000          -0.29000
## 5000 meter buffer          0.46880       -0.1665000          -0.39550
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 1000 meter buffer       -0.4621      0.035280            -0.325200
## 2250 meter buffer       -0.2772      0.668800            -0.114300
## 2500 meter buffer       -0.3628      0.715700            -0.093810
## 750 meter buffer        -0.4509      0.004044            -0.291400
## 500 meter buffer        -0.3458     -0.025160            -0.286300
## 2000 meter buffer       -0.2613      0.587700            -0.147200
## 2750 meter buffer       -0.3731      0.713500            -0.091470
## 3000 meter buffer       -0.3673      0.759500            -0.071460
## 1250 meter buffer       -0.2649      0.138600            -0.311700
## 3250 meter buffer       -0.3980      0.709900            -0.060420
## 4250 meter buffer       -0.7518      0.487200             0.001271
## 4500 meter buffer       -0.7742      0.379300             0.006289
## 1500 meter buffer       -0.2519      0.489900            -0.245300
## 4750 meter buffer       -0.7550      0.385200            -0.015980
## 1750 meter buffer       -0.1555      0.455900            -0.216000
## 4000 meter buffer       -0.7350      0.513400            -0.007683
## 250 meter buffer        -0.0961     -0.311300            -0.199800
## 3500 meter buffer       -0.4239      0.628500            -0.057820
## 3750 meter buffer       -0.5960      0.597500            -0.034360
## 5000 meter buffer       -0.6737      0.207600            -0.063910
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 1000 meter buffer          0.046110       0.11410          -0.12310
## 2250 meter buffer          0.025740      -0.05722          -0.36660
## 2500 meter buffer         -0.023460      -0.05552          -0.26180
## 750 meter buffer           0.026370       0.08005          -0.05458
## 500 meter buffer          -0.048660       0.05504          -0.03073
## 2000 meter buffer          0.014350      -0.05207          -0.34080
## 2750 meter buffer         -0.077060      -0.07645          -0.22520
## 3000 meter buffer         -0.075430      -0.09049          -0.17290
## 1250 meter buffer          0.016590       0.07801          -0.24080
## 3250 meter buffer         -0.055620      -0.09053          -0.07487
## 4250 meter buffer          0.121400      -0.05582           0.21500
## 4500 meter buffer          0.125000      -0.04925           0.22600
## 1500 meter buffer          0.004823       0.06190          -0.21700
## 4750 meter buffer          0.137500      -0.01248           0.23230
## 1750 meter buffer         -0.005642       0.05866          -0.31110
## 4000 meter buffer          0.111700      -0.03095           0.17610
## 250 meter buffer          -0.094330      -0.06221          -0.04710
## 3500 meter buffer         -0.012590      -0.03277          -0.01026
## 3750 meter buffer          0.059180      -0.01418           0.11130
## 5000 meter buffer          0.117500       0.01003           0.23000
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 1000 meter buffer           0.08012      0.034220 18 -268.012 575.2  0.00
## 2250 meter buffer          -0.13470     -0.117200 18 -268.432 576.1  0.84
## 2500 meter buffer          -0.24200     -0.005131 18 -268.638 576.5  1.25
## 750 meter buffer            0.04994      0.071260 18 -268.830 576.9  1.64
## 500 meter buffer            0.06297      0.107700 18 -269.034 577.3  2.04
## 2000 meter buffer          -0.13190     -0.075490 18 -269.328 577.9  2.63
## 2750 meter buffer          -0.38770      0.047000 18 -269.893 579.0  3.76
## 3000 meter buffer          -0.53950      0.052000 18 -270.558 580.3  5.09
## 1250 meter buffer           0.12040     -0.035230 18 -270.684 580.6  5.34
## 3250 meter buffer          -0.61820      0.169200 18 -271.005 581.2  5.99
## 4250 meter buffer          -0.88250      0.595000 18 -271.430 582.1  6.84
## 4500 meter buffer          -0.89610      0.649600 18 -271.748 582.7  7.47
## 1500 meter buffer          -0.11120     -0.086720 18 -271.945 583.1  7.86
## 4750 meter buffer          -0.85780      0.668100 18 -272.202 583.6  8.38
## 1750 meter buffer          -0.02368     -0.084710 18 -272.209 583.6  8.39
## 4000 meter buffer          -0.81260      0.515800 18 -272.249 583.7  8.47
## 250 meter buffer           -0.11080     -0.104400 18 -272.389 584.0  8.75
## 3500 meter buffer          -0.62960      0.231700 18 -272.704 584.6  9.38
## 3750 meter buffer          -0.76020      0.389200 18 -272.890 585.0  9.76
## 5000 meter buffer          -0.79070      0.643000 18 -273.735 586.7 11.44
##                   weight
## 1000 meter buffer  0.266
## 2250 meter buffer  0.175
## 2500 meter buffer  0.142
## 750 meter buffer   0.117
## 500 meter buffer   0.096
## 2000 meter buffer  0.071
## 2750 meter buffer  0.040
## 3000 meter buffer  0.021
## 1250 meter buffer  0.018
## 3250 meter buffer  0.013
## 4250 meter buffer  0.009
## 4500 meter buffer  0.006
## 1500 meter buffer  0.005
## 4750 meter buffer  0.004
## 1750 meter buffer  0.004
## 4000 meter buffer  0.004
## 250 meter buffer   0.003
## 3500 meter buffer  0.002
## 3750 meter buffer  0.002
## 5000 meter buffer  0.001
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For fisher top model is 1000 meter

Model summary 1000m

Let’s print the summary for this model

summary(fisher_mods$`1000 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(fisher, absent_fisher) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    572.0    634.1   -268.0    536.0      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.2788   0.528   
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -3.10303    0.24991 -12.416  < 2e-16 ***
## scale(harvest)            -0.21299    0.11211  -1.900  0.05745 .  
## scale(pipeline)           -0.46212    0.21441  -2.155  0.03114 *  
## scale(roads)               0.03528    0.23091   0.153  0.87858    
## scale(seismic_lines)       0.04611    0.09830   0.469  0.63900    
## scale(seismic_lines_3D)   -0.32519    0.17061  -1.906  0.05665 .  
## scale(trails)              0.11406    0.09359   1.219  0.22296    
## scale(transmission_lines) -0.12312    0.19015  -0.648  0.51730    
## scale(veg_edges)           0.08012    0.21226   0.377  0.70584    
## scale(wells)               0.03422    0.12323   0.278  0.78126    
## scale(osm_industrial)     -0.08007    0.15130  -0.529  0.59666    
## scale(lc_grassland)        0.17267    0.13102   1.318  0.18754    
## scale(lc_coniferous)       0.98548    0.39150   2.517  0.01183 *  
## scale(lc_broadleaf)        0.72202    0.27233   2.651  0.00802 ** 
## scale(lc_mixed)            0.43743    0.24507   1.785  0.07427 .  
## scale(lc_developed)        0.49198    0.17587   2.797  0.00515 ** 
## scale(lc_shrub)           -0.15344    0.26651  -0.576  0.56479    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Everything looks good

Subset models

Add later if deemed necessary.

Top model

Grey wolf

Global models

wolf_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(grey_wolf, absent_grey_wolf) ~
                        
                         # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(wolf_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 500 meter buffer      -3.231          +     -0.037060        -0.971700
## 3500 meter buffer     -3.238          +      0.217800         0.046230
## 3750 meter buffer     -3.238          +      0.235300         0.126700
## 4000 meter buffer     -3.232          +      0.248500         0.224800
## 3250 meter buffer     -3.234          +      0.200700        -0.001385
## 750 meter buffer      -3.245          +      0.037480        -0.889900
## 4250 meter buffer     -3.222          +      0.219000         0.340700
## 3000 meter buffer     -3.230          +      0.188000        -0.074530
## 250 meter buffer      -3.268          +     -0.027310        -0.672400
## 2000 meter buffer     -3.252          +     -0.008536        -0.424000
## 1750 meter buffer     -3.250          +     -0.018120        -0.390000
## 2750 meter buffer     -3.229          +      0.156400        -0.142000
## 1000 meter buffer     -3.232          +      0.024070        -0.524800
## 1500 meter buffer     -3.240          +     -0.012710        -0.334600
## 2250 meter buffer     -3.241          +      0.028360        -0.380400
## 4500 meter buffer     -3.210          +      0.179600         0.387400
## 2500 meter buffer     -3.228          +      0.093950        -0.240700
## 4750 meter buffer     -3.208          +      0.149400         0.408200
## 5000 meter buffer     -3.205          +      0.166600         0.400900
## 1250 meter buffer     -3.219          +      0.027900        -0.350400
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 500 meter buffer          -0.98650        -0.048790        -0.313100
## 3500 meter buffer          0.39970         0.165800         0.580500
## 3750 meter buffer          0.48230         0.169400         0.607300
## 4000 meter buffer          0.60450         0.150700         0.629700
## 3250 meter buffer          0.34270         0.198500         0.499300
## 750 meter buffer          -0.89080        -0.007924        -0.149800
## 4250 meter buffer          0.75150         0.156800         0.608400
## 3000 meter buffer          0.25060         0.214200         0.420900
## 250 meter buffer          -0.51490        -0.233900        -0.221000
## 2000 meter buffer         -0.10830         0.328200         0.090040
## 1750 meter buffer         -0.09863         0.349900         0.050860
## 2750 meter buffer          0.19760         0.256000         0.341900
## 1000 meter buffer         -0.40450         0.242500        -0.062530
## 1500 meter buffer         -0.06133         0.364000         0.013020
## 2250 meter buffer         -0.05733         0.311200         0.138400
## 4500 meter buffer          0.80090         0.083230         0.589100
## 2500 meter buffer          0.09298         0.314300         0.227500
## 4750 meter buffer          0.82190         0.028690         0.589600
## 5000 meter buffer          0.80830         0.062450         0.598900
## 1250 meter buffer         -0.14750         0.294400         0.005741
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 500 meter buffer        -0.4470000         -0.48920         -0.047730
## 3500 meter buffer       -0.0122800          0.12360         -0.137300
## 3750 meter buffer       -0.0082400          0.18980         -0.102700
## 4000 meter buffer        0.0417800          0.27260         -0.025010
## 3250 meter buffer       -0.0100600          0.09504         -0.117900
## 750 meter buffer        -0.3254000         -0.51950          0.008087
## 4250 meter buffer        0.1206000          0.37130          0.011010
## 3000 meter buffer       -0.0133700          0.07169         -0.091320
## 250 meter buffer        -0.1315000         -0.31460          0.083260
## 2000 meter buffer        0.0127700         -0.27740         -0.187200
## 1750 meter buffer        0.0216700         -0.29160         -0.136500
## 2750 meter buffer        0.0069020          0.03526         -0.086450
## 1000 meter buffer       -0.1004000         -0.29860         -0.002024
## 1500 meter buffer        0.0376400         -0.25910         -0.097470
## 2250 meter buffer       -0.0097320         -0.21460         -0.202000
## 4500 meter buffer        0.1565000          0.38040         -0.052350
## 2500 meter buffer       -0.0003944         -0.05267         -0.143300
## 4750 meter buffer        0.1515000          0.33260         -0.201000
## 5000 meter buffer        0.1053000          0.30490         -0.230100
## 1250 meter buffer        0.0059370         -0.18650         -0.019950
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 500 meter buffer        0.34890      -0.38150              -0.9052
## 3500 meter buffer      -0.44600      -0.02598              -0.7740
## 3750 meter buffer      -0.46090      -0.05365              -0.7688
## 4000 meter buffer      -0.52120      -0.12340              -0.7648
## 3250 meter buffer      -0.35390      -0.03936              -0.8080
## 750 meter buffer        0.38540      -0.58220              -0.9829
## 4250 meter buffer      -0.43160      -0.12370              -0.7782
## 3000 meter buffer      -0.31250      -0.08431              -0.8640
## 250 meter buffer        0.22500      -0.18480              -1.2650
## 2000 meter buffer       0.19270      -0.33400              -0.9736
## 1750 meter buffer       0.29720      -0.46840              -0.9661
## 2750 meter buffer      -0.24670      -0.10280              -0.9034
## 1000 meter buffer       0.37390      -0.67980              -0.9062
## 1500 meter buffer       0.34370      -0.62370              -0.9327
## 2250 meter buffer       0.03818      -0.14240              -0.9378
## 4500 meter buffer      -0.35220       0.01628              -0.7679
## 2500 meter buffer      -0.11820      -0.10630              -0.9231
## 4750 meter buffer      -0.40130       0.25730              -0.6983
## 5000 meter buffer      -0.52980       0.31470              -0.6436
## 1250 meter buffer       0.33120      -0.63520              -0.9027
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 500 meter buffer           0.125000       0.08796          -0.19730
## 3500 meter buffer         -0.089180       0.06981           0.26460
## 3750 meter buffer         -0.100500       0.06518           0.26380
## 4000 meter buffer         -0.103300       0.03829           0.29540
## 3250 meter buffer         -0.063300       0.07017           0.22960
## 750 meter buffer           0.157100      -0.01028          -0.12220
## 4250 meter buffer         -0.113600       0.05983           0.20550
## 3000 meter buffer         -0.029190       0.05093           0.24090
## 250 meter buffer           0.036700       0.20850          -0.25300
## 2000 meter buffer          0.021150       0.21190          -0.05804
## 1750 meter buffer         -0.018600       0.21330          -0.11380
## 2750 meter buffer         -0.003554       0.05788           0.21810
## 1000 meter buffer          0.086680       0.07459          -0.06863
## 1500 meter buffer         -0.002519       0.20930          -0.13230
## 2250 meter buffer          0.024590       0.16410           0.02337
## 4500 meter buffer         -0.130500       0.08058           0.08265
## 2500 meter buffer          0.015410       0.10940           0.14540
## 4750 meter buffer         -0.131200       0.07740           0.00401
## 5000 meter buffer         -0.114900       0.05334           0.05140
## 1250 meter buffer          0.024350       0.13480          -0.09049
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 500 meter buffer          -0.145700      -0.28210 18 -236.653 512.5  0.00
## 3500 meter buffer         -0.188000      -0.22190 18 -237.226 513.7  1.15
## 3750 meter buffer         -0.176000      -0.24700 18 -237.463 514.1  1.62
## 4000 meter buffer         -0.176000      -0.16900 18 -237.738 514.7  2.17
## 3250 meter buffer         -0.236700      -0.21760 18 -238.268 515.7  3.23
## 750 meter buffer          -0.092290      -0.02706 18 -238.652 516.5  4.00
## 4250 meter buffer         -0.156000      -0.21800 18 -239.238 517.7  5.17
## 3000 meter buffer         -0.284700      -0.14220 18 -239.529 518.3  5.75
## 250 meter buffer          -0.082900      -0.06612 18 -239.769 518.7  6.23
## 2000 meter buffer         -0.156100      -0.13170 18 -239.864 518.9  6.42
## 1750 meter buffer         -0.093230      -0.12370 18 -239.948 519.1  6.59
## 2750 meter buffer         -0.331100      -0.09790 18 -240.224 519.7  7.14
## 1000 meter buffer         -0.124800      -0.06659 18 -240.303 519.8  7.30
## 1500 meter buffer         -0.032900      -0.06583 18 -240.428 520.1  7.55
## 2250 meter buffer         -0.240100      -0.19080 18 -240.681 520.6  8.06
## 4500 meter buffer         -0.085510      -0.31020 18 -240.762 520.7  8.22
## 2500 meter buffer         -0.323900      -0.14500 18 -241.169 521.5  9.03
## 4750 meter buffer          0.005409      -0.44860 18 -241.369 521.9  9.43
## 5000 meter buffer         -0.002138      -0.45350 18 -241.521 522.3  9.74
## 1250 meter buffer         -0.102500      -0.03946 18 -241.850 522.9 10.39
##                   weight
## 500 meter buffer   0.325
## 3500 meter buffer  0.183
## 3750 meter buffer  0.145
## 4000 meter buffer  0.110
## 3250 meter buffer  0.065
## 750 meter buffer   0.044
## 4250 meter buffer  0.024
## 3000 meter buffer  0.018
## 250 meter buffer   0.014
## 2000 meter buffer  0.013
## 1750 meter buffer  0.012
## 2750 meter buffer  0.009
## 1000 meter buffer  0.008
## 1500 meter buffer  0.007
## 2250 meter buffer  0.006
## 4500 meter buffer  0.005
## 2500 meter buffer  0.004
## 4750 meter buffer  0.003
## 5000 meter buffer  0.002
## 1250 meter buffer  0.002
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For grey wolf top model is 500 m buffer

Model summary 500m

Let’s get the model summary for this model

summary(wolf_mods$`500 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(grey_wolf, absent_grey_wolf) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    509.3    571.3   -236.7    473.3      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance  Std.Dev. 
##  array  (Intercept) 5.924e-10 2.434e-05
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -3.23077    0.13696 -23.589  < 2e-16 ***
## scale(harvest)            -0.03706    0.09169  -0.404  0.68611    
## scale(pipeline)            0.34887    0.15217   2.293  0.02187 *  
## scale(roads)              -0.38147    0.19999  -1.907  0.05646 .  
## scale(seismic_lines)       0.12504    0.08754   1.428  0.15320    
## scale(seismic_lines_3D)   -0.90518    0.27643  -3.275  0.00106 ** 
## scale(trails)              0.08796    0.08212   1.071  0.28410    
## scale(transmission_lines) -0.19728    0.14199  -1.389  0.16470    
## scale(veg_edges)          -0.14566    0.14749  -0.988  0.32336    
## scale(wells)              -0.28215    0.15563  -1.813  0.06984 .  
## scale(osm_industrial)     -0.04773    0.10158  -0.470  0.63844    
## scale(lc_grassland)       -0.31308    0.15267  -2.051  0.04030 *  
## scale(lc_coniferous)      -0.98649    0.41080  -2.401  0.01633 *  
## scale(lc_broadleaf)       -0.97169    0.29859  -3.254  0.00114 ** 
## scale(lc_mixed)           -0.44697    0.24027  -1.860  0.06284 .  
## scale(lc_developed)       -0.04879    0.17199  -0.284  0.77665    
## scale(lc_shrub)           -0.48916    0.23223  -2.106  0.03517 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

interesting that seismic lines weren’t significant but 3d seismic were and have a negative estimate

Subset models

Add later if deemed necessary.

Top model

Lynx

Global models

lynx_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(lynx, absent_lynx) ~
                        
                       # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(lynx_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 1250 meter buffer     -1.998          +      -0.10700         0.951500
## 1000 meter buffer     -1.999          +      -0.10950         1.200000
## 1500 meter buffer     -1.983          +      -0.08890         0.510000
## 500 meter buffer      -1.979          +      -0.20620         1.094000
## 750 meter buffer      -1.982          +      -0.14150         1.202000
## 1750 meter buffer     -1.972          +      -0.04953         0.395600
## 250 meter buffer      -1.957          +      -0.31350         0.514700
## 2000 meter buffer     -1.964          +      -0.01751         0.297200
## 2250 meter buffer     -1.951          +       0.03254         0.173700
## 2500 meter buffer     -1.940          +       0.08284         0.065000
## 2750 meter buffer     -1.934          +       0.12560        -0.002596
## 5000 meter buffer     -1.920          +       0.15490        -0.271800
## 3000 meter buffer     -1.932          +       0.14720        -0.008921
## 3250 meter buffer     -1.931          +       0.15490        -0.009491
## 4750 meter buffer     -1.919          +       0.16900        -0.306100
## 3500 meter buffer     -1.926          +       0.15950        -0.064560
## 3750 meter buffer     -1.925          +       0.16140        -0.139600
## 4250 meter buffer     -1.923          +       0.17840        -0.249800
## 4000 meter buffer     -1.923          +       0.16860        -0.220800
## 4500 meter buffer     -1.920          +       0.17420        -0.271800
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 1250 meter buffer           2.2940           0.5731          0.52720
## 1000 meter buffer           2.5530           0.6433          0.59470
## 1500 meter buffer           1.6460           0.4934          0.36640
## 500 meter buffer            2.2360           0.5143          0.72110
## 750 meter buffer            2.4650           0.6015          0.62100
## 1750 meter buffer           1.4670           0.4982          0.26500
## 250 meter buffer            1.3370           0.2788          0.44470
## 2000 meter buffer           1.3040           0.4947          0.23210
## 2250 meter buffer           1.0990           0.4593          0.18130
## 2500 meter buffer           0.8984           0.4509          0.16350
## 2750 meter buffer           0.7380           0.4410          0.18500
## 5000 meter buffer           0.3504           0.3159         -0.02850
## 3000 meter buffer           0.6937           0.4508          0.18720
## 3250 meter buffer           0.7028           0.4774          0.19180
## 4750 meter buffer           0.3355           0.3516         -0.03058
## 3500 meter buffer           0.6110           0.4454          0.15500
## 3750 meter buffer           0.5303           0.4425          0.11400
## 4250 meter buffer           0.4247           0.4277          0.06568
## 4000 meter buffer           0.4381           0.4143          0.07779
## 4500 meter buffer           0.3901           0.4032          0.01749
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 1250 meter buffer          1.26600           1.2370           0.18780
## 1000 meter buffer          1.41900           1.3690           0.22610
## 1500 meter buffer          0.89640           0.8896           0.12510
## 500 meter buffer           1.14200           1.1440           0.22890
## 750 meter buffer           1.35700           1.2830           0.25180
## 1750 meter buffer          0.77230           0.7947           0.10610
## 250 meter buffer           0.73140           0.7349           0.16530
## 2000 meter buffer          0.67270           0.6892           0.08232
## 2250 meter buffer          0.53040           0.5593           0.06887
## 2500 meter buffer          0.39370           0.4528           0.07333
## 2750 meter buffer          0.26820           0.3755           0.05677
## 5000 meter buffer          0.09991           0.1316          -0.05720
## 3000 meter buffer          0.21590           0.3483           0.05043
## 3250 meter buffer          0.22290           0.3610           0.06459
## 4750 meter buffer          0.09305           0.1285          -0.02188
## 3500 meter buffer          0.17910           0.3089           0.03058
## 3750 meter buffer          0.14790           0.2674           0.02198
## 4250 meter buffer          0.12320           0.2027           0.03087
## 4000 meter buffer          0.11450           0.2126           0.01585
## 4500 meter buffer          0.11010           0.1744           0.02351
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 1250 meter buffer      0.011650       0.18440              0.08558
## 1000 meter buffer      0.044020       0.04073              0.08290
## 1500 meter buffer      0.018580       0.34000              0.10820
## 500 meter buffer      -0.008030       0.03130              0.06549
## 750 meter buffer      -0.003399       0.01784              0.06542
## 1750 meter buffer      0.014920       0.46950              0.12480
## 250 meter buffer       0.016210       0.02852              0.04631
## 2000 meter buffer     -0.014960       0.53310              0.13510
## 2250 meter buffer     -0.081600       0.56500              0.13330
## 2500 meter buffer     -0.107300       0.42190              0.12300
## 2750 meter buffer     -0.077360       0.34500              0.12380
## 5000 meter buffer      0.082550       0.39330              0.14820
## 3000 meter buffer     -0.061900       0.29280              0.12390
## 3250 meter buffer     -0.044800       0.21130              0.11950
## 4750 meter buffer      0.022940       0.30420              0.11610
## 3500 meter buffer     -0.061890       0.23630              0.10620
## 3750 meter buffer     -0.039720       0.20790              0.09828
## 4250 meter buffer     -0.028600       0.15290              0.09666
## 4000 meter buffer     -0.033380       0.18580              0.09580
## 4500 meter buffer      0.008188       0.19140              0.10340
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 1250 meter buffer         -0.003312       0.11570         -0.083150
## 1000 meter buffer         -0.031170       0.10750         -0.047370
## 1500 meter buffer         -0.007120       0.12780         -0.075910
## 500 meter buffer          -0.086320      -0.02598         -0.078670
## 750 meter buffer          -0.044540       0.02683         -0.004659
## 1750 meter buffer         -0.028230       0.11330         -0.047020
## 250 meter buffer          -0.153400       0.03605         -0.028840
## 2000 meter buffer         -0.026320       0.11000         -0.023570
## 2250 meter buffer         -0.055640       0.09932          0.039340
## 2500 meter buffer         -0.072270       0.08368          0.087040
## 2750 meter buffer         -0.075140       0.08633          0.053090
## 5000 meter buffer         -0.123300       0.13080         -0.097910
## 3000 meter buffer         -0.089490       0.09532          0.021930
## 3250 meter buffer         -0.109400       0.09882          0.012100
## 4750 meter buffer         -0.115700       0.12240         -0.025950
## 3500 meter buffer         -0.123300       0.10480          0.035380
## 3750 meter buffer         -0.126600       0.10650          0.023980
## 4250 meter buffer         -0.123400       0.10910          0.026350
## 4000 meter buffer         -0.124500       0.10590          0.018090
## 4500 meter buffer         -0.110100       0.11820         -0.009629
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 1250 meter buffer           -0.3415       -0.2457 18 -414.250 867.7  0.00
## 1000 meter buffer           -0.3337       -0.1743 18 -414.822 868.9  1.14
## 1500 meter buffer           -0.4351       -0.2975 18 -416.528 872.3  4.56
## 500 meter buffer            -0.3271       -0.1569 18 -416.655 872.5  4.81
## 750 meter buffer            -0.3294       -0.1261 18 -418.644 876.5  8.79
## 1750 meter buffer           -0.5412       -0.3438 18 -419.271 877.8 10.04
## 250 meter buffer            -0.1747       -0.1823 18 -420.007 879.2 11.51
## 2000 meter buffer           -0.5642       -0.4024 18 -421.235 881.7 13.97
## 2250 meter buffer           -0.5536       -0.4200 18 -424.464 888.1 20.43
## 2500 meter buffer           -0.4965       -0.3664 18 -427.631 894.5 26.76
## 2750 meter buffer           -0.4432       -0.3658 18 -429.440 898.1 30.38
## 5000 meter buffer           -0.2836       -0.3014 18 -429.999 899.2 31.50
## 3000 meter buffer           -0.4081       -0.3515 18 -430.196 899.6 31.89
## 3250 meter buffer           -0.3866       -0.3170 18 -430.805 900.8 33.11
## 4750 meter buffer           -0.3000       -0.2467 18 -431.395 902.0 34.29
## 3500 meter buffer           -0.3665       -0.2986 18 -431.591 902.4 34.68
## 3750 meter buffer           -0.3423       -0.2757 18 -432.198 903.6 35.90
## 4250 meter buffer           -0.3173       -0.2238 18 -432.979 905.2 37.46
## 4000 meter buffer           -0.3108       -0.2373 18 -433.147 905.5 37.79
## 4500 meter buffer           -0.3208       -0.2096 18 -433.179 905.6 37.86
##                   weight
## 1250 meter buffer  0.562
## 1000 meter buffer  0.317
## 1500 meter buffer  0.058
## 500 meter buffer   0.051
## 750 meter buffer   0.007
## 1750 meter buffer  0.004
## 250 meter buffer   0.002
## 2000 meter buffer  0.001
## 2250 meter buffer  0.000
## 2500 meter buffer  0.000
## 2750 meter buffer  0.000
## 5000 meter buffer  0.000
## 3000 meter buffer  0.000
## 3250 meter buffer  0.000
## 4750 meter buffer  0.000
## 3500 meter buffer  0.000
## 3750 meter buffer  0.000
## 4250 meter buffer  0.000
## 4000 meter buffer  0.000
## 4500 meter buffer  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For lynx the top model is the 1250 m buffer

Model summary 1250m

Let’s get the model summary

summary(lynx_mods$`1250 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(lynx, absent_lynx) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    864.5    926.5   -414.2    828.5      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.02967  0.1722  
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.997708   0.097494 -20.491  < 2e-16 ***
## scale(harvest)            -0.106953   0.081260  -1.316  0.18812    
## scale(pipeline)            0.011650   0.114082   0.102  0.91866    
## scale(roads)               0.184363   0.187036   0.986  0.32428    
## scale(seismic_lines)      -0.003312   0.067892  -0.049  0.96110    
## scale(seismic_lines_3D)    0.085578   0.061810   1.385  0.16620    
## scale(trails)              0.115665   0.055589   2.081  0.03746 *  
## scale(transmission_lines) -0.083153   0.106651  -0.780  0.43558    
## scale(veg_edges)          -0.341523   0.134226  -2.544  0.01095 *  
## scale(wells)              -0.245720   0.098236  -2.501  0.01237 *  
## scale(osm_industrial)      0.187831   0.072422   2.594  0.00950 ** 
## scale(lc_grassland)        0.527181   0.127159   4.146 3.39e-05 ***
## scale(lc_coniferous)       2.294419   0.554694   4.136 3.53e-05 ***
## scale(lc_broadleaf)        0.951507   0.358493   2.654  0.00795 ** 
## scale(lc_mixed)            1.265559   0.314719   4.021 5.79e-05 ***
## scale(lc_developed)        0.573095   0.144770   3.959 7.54e-05 ***
## scale(lc_shrub)            1.237013   0.306204   4.040 5.35e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Everything looks good so far

Subset models

Add later if deemed necessary.

Top model

Moose

Global model

moose_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(moose, absent_moose) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))

# model selection
model.sel(moose_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 250 meter buffer      -1.878          +     -0.002032         -0.04093
## 500 meter buffer      -1.881          +      0.007602         -0.15570
## 750 meter buffer      -1.879          +      0.018560          0.04249
## 5000 meter buffer     -1.882          +      0.005041          0.97770
## 1000 meter buffer     -1.871          +      0.010350          0.00255
## 4750 meter buffer     -1.877          +      0.002859          0.93390
## 3750 meter buffer     -1.869          +      0.023660          0.60380
## 4000 meter buffer     -1.869          +      0.022960          0.69640
## 4500 meter buffer     -1.869          +      0.002659          0.85200
## 3500 meter buffer     -1.865          +      0.033830          0.54850
## 4250 meter buffer     -1.867          +      0.010490          0.77670
## 3250 meter buffer     -1.862          +      0.046550          0.47340
## 3000 meter buffer     -1.863          +      0.057610          0.40840
## 2750 meter buffer     -1.860          +      0.064590          0.32600
## 2500 meter buffer     -1.856          +      0.070920          0.28420
## 1250 meter buffer     -1.857          +      0.039370         -0.01379
## 2250 meter buffer     -1.849          +      0.067100          0.21110
## 1500 meter buffer     -1.851          +      0.042190         -0.01124
## 1750 meter buffer     -1.845          +      0.041230          0.03104
## 2000 meter buffer     -1.843          +      0.032260          0.13520
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 250 meter buffer          -0.58030         -0.33560       -0.1402000
## 500 meter buffer          -0.74630         -0.18660       -0.1253000
## 750 meter buffer          -0.47410         -0.08463       -0.0381800
## 5000 meter buffer          0.79260          0.54480        0.0062190
## 1000 meter buffer         -0.50940         -0.08745       -0.0597700
## 4750 meter buffer          0.72150          0.51900        0.0475100
## 3750 meter buffer          0.19970          0.24010       -0.0047560
## 4000 meter buffer          0.32480          0.35190       -0.0001991
## 4500 meter buffer          0.57510          0.41390        0.0439300
## 3500 meter buffer          0.12590          0.18550        0.0138500
## 4250 meter buffer          0.45770          0.37710        0.0269800
## 3250 meter buffer          0.03208          0.09211        0.0338300
## 3000 meter buffer         -0.05974         -0.01656        0.0534800
## 2750 meter buffer         -0.18380         -0.14370        0.0713600
## 2500 meter buffer         -0.23100         -0.13840        0.0705900
## 1250 meter buffer         -0.52390         -0.07189       -0.0273000
## 2250 meter buffer         -0.31020         -0.10800        0.0484400
## 1500 meter buffer         -0.51040         -0.04036       -0.0108500
## 1750 meter buffer         -0.47420         -0.09248        0.0164000
## 2000 meter buffer         -0.36670         -0.10310        0.0273400
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 250 meter buffer          -0.27790        -0.088590           -0.1908
## 500 meter buffer          -0.38860        -0.236700           -0.3299
## 750 meter buffer          -0.17030        -0.134300           -0.3327
## 5000 meter buffer          0.21320         0.604400           -0.2861
## 1000 meter buffer         -0.20600        -0.126800           -0.2908
## 4750 meter buffer          0.20820         0.563100           -0.2792
## 3750 meter buffer          0.08199         0.286600           -0.4126
## 4000 meter buffer          0.11870         0.365200           -0.3352
## 4500 meter buffer          0.17090         0.489900           -0.3017
## 3500 meter buffer          0.04727         0.232000           -0.3975
## 4250 meter buffer          0.15780         0.431400           -0.3017
## 3250 meter buffer         -0.01014         0.157900           -0.4158
## 3000 meter buffer         -0.07013         0.099760           -0.4314
## 2750 meter buffer         -0.15070         0.029640           -0.4406
## 2500 meter buffer         -0.17710         0.009509           -0.3517
## 1250 meter buffer         -0.24480        -0.125200           -0.2146
## 2250 meter buffer         -0.19270        -0.026350           -0.2546
## 1500 meter buffer         -0.26390        -0.135500           -0.2362
## 1750 meter buffer         -0.26020        -0.125100           -0.2128
## 2000 meter buffer         -0.19500        -0.076940           -0.2405
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 250 meter buffer      -0.006781      -0.19920              0.18590
## 500 meter buffer       0.102800      -0.22250              0.15950
## 750 meter buffer       0.174600      -0.22750              0.12520
## 5000 meter buffer     -0.127700       0.33710              0.05034
## 1000 meter buffer      0.232800      -0.25860              0.08750
## 4750 meter buffer     -0.132200       0.31090              0.06682
## 3750 meter buffer     -0.010510       0.38230              0.03537
## 4000 meter buffer     -0.003876       0.27980              0.04015
## 4500 meter buffer     -0.083720       0.30210              0.05842
## 3500 meter buffer     -0.027980       0.37520              0.03467
## 4250 meter buffer     -0.066350       0.25880              0.04796
## 3250 meter buffer     -0.065050       0.43670              0.03313
## 3000 meter buffer     -0.101600       0.49180              0.02630
## 2750 meter buffer     -0.131500       0.49970              0.01676
## 2500 meter buffer     -0.160600       0.34010              0.01979
## 1250 meter buffer      0.186900      -0.28900              0.07527
## 2250 meter buffer     -0.099920       0.13580              0.03606
## 1500 meter buffer      0.132800      -0.25840              0.07116
## 1750 meter buffer      0.040180      -0.11020              0.05845
## 2000 meter buffer      0.006445       0.03398              0.04289
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 250 meter buffer           0.041250      0.052820         9.152e-02
## 500 meter buffer           0.028140      0.045090         3.699e-02
## 750 meter buffer          -0.025640      0.035850        -8.078e-03
## 5000 meter buffer         -0.078360     -0.015410         3.037e-01
## 1000 meter buffer         -0.077310      0.045370        -5.423e-02
## 4750 meter buffer         -0.081570     -0.008001         2.681e-01
## 3750 meter buffer         -0.100500      0.037600         2.074e-01
## 4000 meter buffer         -0.107100      0.025630         2.002e-01
## 4500 meter buffer         -0.091260      0.010430         2.228e-01
## 3500 meter buffer         -0.097610      0.042530         2.104e-01
## 4250 meter buffer         -0.104300      0.013880         2.196e-01
## 3250 meter buffer         -0.074480      0.044980         2.194e-01
## 3000 meter buffer         -0.057350      0.047050         2.386e-01
## 2750 meter buffer         -0.022730      0.051450         2.552e-01
## 2500 meter buffer         -0.009657      0.044060         2.791e-01
## 1250 meter buffer         -0.036630      0.047520        -2.164e-02
## 2250 meter buffer         -0.013640      0.030010         2.044e-01
## 1500 meter buffer         -0.042620      0.061730         6.225e-05
## 1750 meter buffer         -0.044150      0.043640         7.419e-02
## 2000 meter buffer         -0.019850      0.049330         8.495e-02
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 250 meter buffer           0.110900      -0.05671 18 -408.795 856.8  0.00
## 500 meter buffer           0.011790      -0.07600 18 -408.932 857.1  0.27
## 750 meter buffer          -0.010410      -0.09323 18 -411.103 861.4  4.62
## 5000 meter buffer         -0.626700      -0.22700 18 -413.830 866.9 10.07
## 1000 meter buffer          0.017430      -0.09916 18 -414.280 867.8 10.97
## 4750 meter buffer         -0.575100      -0.26440 18 -416.094 871.4 14.60
## 3750 meter buffer         -0.390600      -0.27620 18 -417.695 874.6 17.80
## 4000 meter buffer         -0.460500      -0.24140 18 -417.747 874.7 17.90
## 4500 meter buffer         -0.498600      -0.24690 18 -417.903 875.0 18.22
## 3500 meter buffer         -0.343500      -0.30270 18 -418.307 875.8 19.02
## 4250 meter buffer         -0.463700      -0.21590 18 -418.360 875.9 19.13
## 3250 meter buffer         -0.282800      -0.33550 18 -419.112 877.4 20.63
## 3000 meter buffer         -0.212500      -0.36570 18 -419.159 877.5 20.73
## 2750 meter buffer         -0.115200      -0.34760 18 -419.237 877.7 20.88
## 2500 meter buffer         -0.056660      -0.32340 18 -419.865 878.9 22.14
## 1250 meter buffer          0.042460      -0.06342 18 -419.892 879.0 22.19
## 2250 meter buffer          0.006556      -0.28980 18 -422.400 884.0 27.21
## 1500 meter buffer          0.085920      -0.11490 18 -422.726 884.7 27.86
## 1750 meter buffer          0.060070      -0.20410 18 -424.194 887.6 30.80
## 2000 meter buffer          0.082470      -0.29880 18 -424.559 888.3 31.53
##                   weight
## 250 meter buffer   0.504
## 500 meter buffer   0.440
## 750 meter buffer   0.050
## 5000 meter buffer  0.003
## 1000 meter buffer  0.002
## 4750 meter buffer  0.000
## 3750 meter buffer  0.000
## 4000 meter buffer  0.000
## 4500 meter buffer  0.000
## 3500 meter buffer  0.000
## 4250 meter buffer  0.000
## 3250 meter buffer  0.000
## 3000 meter buffer  0.000
## 2750 meter buffer  0.000
## 2500 meter buffer  0.000
## 1250 meter buffer  0.000
## 2250 meter buffer  0.000
## 1500 meter buffer  0.000
## 1750 meter buffer  0.000
## 2000 meter buffer  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For moose the top model is the 250 m buffer

Model summary 250m

Let’s get the model summary

summary(moose_mods$`250 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(moose, absent_moose) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    853.6    915.6   -408.8    817.6      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 0.3108   0.5575  
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -1.877557   0.236425  -7.941    2e-15 ***
## scale(harvest)            -0.002032   0.062622  -0.032 0.974117    
## scale(pipeline)           -0.006781   0.082485  -0.082 0.934477    
## scale(roads)              -0.199217   0.116670  -1.708 0.087723 .  
## scale(seismic_lines)       0.041253   0.063722   0.647 0.517373    
## scale(seismic_lines_3D)    0.185895   0.055703   3.337 0.000846 ***
## scale(trails)              0.052822   0.053881   0.980 0.326916    
## scale(transmission_lines)  0.091525   0.076611   1.195 0.232217    
## scale(veg_edges)           0.110934   0.099668   1.113 0.265694    
## scale(wells)              -0.056712   0.070758  -0.801 0.422849    
## scale(osm_industrial)     -0.190844   0.072207  -2.643 0.008217 ** 
## scale(lc_grassland)       -0.140190   0.106892  -1.312 0.189684    
## scale(lc_coniferous)      -0.580280   0.272816  -2.127 0.033420 *  
## scale(lc_broadleaf)       -0.040927   0.190428  -0.215 0.829830    
## scale(lc_mixed)           -0.277906   0.169142  -1.643 0.100376    
## scale(lc_developed)       -0.335609   0.120540  -2.784 0.005366 ** 
## scale(lc_shrub)           -0.088587   0.152663  -0.580 0.561727    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looks good for now

Subset models

Add later if deemed necessary

Top model

Red fox

Global models

fox_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

     glmmTMB::glmmTMB(cbind(red_fox, absent_red_fox) ~
                        
                         # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))


# model selection
model.sel(fox_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 4750 meter buffer     -5.482          +       0.33900         -0.59670
## 5000 meter buffer     -5.317          +       0.30480         -0.51060
## 4500 meter buffer     -5.305          +       0.34250         -0.83470
## 1000 meter buffer     -4.590          +      -0.23040         -0.27900
## 1250 meter buffer     -4.625          +      -0.25760         -0.29670
## 1500 meter buffer     -4.798          +      -0.18940         -0.15750
## 1750 meter buffer     -4.856          +      -0.11410         -0.04801
## 2500 meter buffer     -4.819          +       0.13900          0.10620
## 2250 meter buffer     -4.687          +      -0.02548          0.11060
## 4250 meter buffer     -4.932          +       0.31360         -1.02500
## 2000 meter buffer     -4.720          +      -0.03125          0.07487
## 750 meter buffer      -4.585          +      -0.23310         -0.10160
## 4000 meter buffer     -4.784          +       0.26090         -0.82660
## 2750 meter buffer     -4.660          +       0.15870         -0.02933
## 3750 meter buffer     -4.717          +       0.15350         -0.65610
## 500 meter buffer      -4.471          +      -0.14050          0.62930
## 3000 meter buffer     -4.645          +       0.18320         -0.24890
## 250 meter buffer      -4.526          +      -0.12290         -0.30710
## 3500 meter buffer     -4.575          +       0.12620         -0.57170
## 3250 meter buffer     -4.501          +       0.22740         -0.63270
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 4750 meter buffer          -1.7340         -1.94600          0.45870
## 5000 meter buffer          -1.6740         -2.14800          0.39790
## 4500 meter buffer          -1.8170         -1.53100          0.42530
## 1000 meter buffer          -0.9658         -0.19720          0.31940
## 1250 meter buffer          -1.0170         -0.30950          0.35830
## 1500 meter buffer          -0.8697         -0.32920          0.52600
## 1750 meter buffer          -0.7501         -0.46890          0.64170
## 2500 meter buffer          -0.6056         -0.70630          0.87200
## 2250 meter buffer          -0.5304         -0.33220          0.87670
## 4250 meter buffer          -1.9290         -1.32700          0.28900
## 2000 meter buffer          -0.5677         -0.40830          0.72750
## 750 meter buffer           -0.7481         -0.01489          0.26660
## 4000 meter buffer          -1.5360         -1.28900          0.31630
## 2750 meter buffer          -0.6417         -0.88420          0.65290
## 3750 meter buffer          -1.2830         -1.28300          0.43740
## 500 meter buffer            0.2211          0.36930          0.26320
## 3000 meter buffer          -0.8074         -0.95970          0.57610
## 250 meter buffer           -1.0680          0.14490         -0.06256
## 3500 meter buffer          -1.1860         -0.91950          0.48310
## 3250 meter buffer          -1.2270         -0.76390          0.52800
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 4750 meter buffer          -2.1350          -0.2193           0.06849
## 5000 meter buffer          -2.0760          -0.1973           0.11500
## 4500 meter buffer          -2.0420          -0.3040           0.05564
## 1000 meter buffer          -0.4364          -0.5876          -0.38600
## 1250 meter buffer          -0.5442          -0.6927          -0.37030
## 1500 meter buffer          -0.5289          -0.6217          -0.22830
## 1750 meter buffer          -0.5698          -0.5523          -0.15290
## 2500 meter buffer          -0.8063          -0.6174          -0.54210
## 2250 meter buffer          -0.5755          -0.6210          -0.50300
## 4250 meter buffer          -1.8380          -0.4017          -0.30790
## 2000 meter buffer          -0.5087          -0.5390          -0.16830
## 750 meter buffer           -0.2679          -0.6480          -0.14420
## 4000 meter buffer          -1.5690          -0.4128          -1.27300
## 2750 meter buffer          -0.8625          -0.4557          -0.33110
## 3750 meter buffer          -1.3670          -0.4864          -1.30000
## 500 meter buffer            0.2701          -0.2161           0.33850
## 3000 meter buffer          -1.0320          -0.5721          -0.46830
## 250 meter buffer           -0.3924          -1.0100           0.08637
## 3500 meter buffer          -1.1970          -0.5730          -0.55310
## 3250 meter buffer          -1.1790          -0.7950          -0.26230
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 4750 meter buffer       0.55450        1.7720             -0.25980
## 5000 meter buffer       0.57000        1.5980             -0.33860
## 4500 meter buffer       0.35170        1.4700             -0.30540
## 1000 meter buffer      -0.08652        0.9211             -0.43730
## 1250 meter buffer       0.13300        0.8973             -0.45770
## 1500 meter buffer      -0.05386        0.6397             -0.36090
## 1750 meter buffer      -0.23330        0.9070             -0.32960
## 2500 meter buffer      -0.66470        1.5980             -0.13210
## 2250 meter buffer      -0.58010        1.2340             -0.26340
## 4250 meter buffer       0.26040        1.2310             -0.49290
## 2000 meter buffer      -0.47410        0.9593             -0.31730
## 750 meter buffer        0.28020        0.5181             -0.25280
## 4000 meter buffer       0.12790        1.8010             -0.41760
## 2750 meter buffer      -0.43360        1.4350             -0.10180
## 3750 meter buffer       0.14840        1.7190             -0.33560
## 500 meter buffer        0.04199        0.1326             -0.17750
## 3000 meter buffer      -0.25270        1.5540             -0.06029
## 250 meter buffer       -0.12160        0.2010              0.02755
## 3500 meter buffer       0.02500        1.1300             -0.14100
## 3250 meter buffer      -0.19720        0.7720             -0.13030
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 4750 meter buffer         -0.778800       -1.1220           0.51430
## 5000 meter buffer         -0.668900       -0.9373           0.37600
## 4500 meter buffer         -0.689600       -1.1380           0.65820
## 1000 meter buffer         -0.034770       -0.6892          -0.84680
## 1250 meter buffer         -0.001617       -0.7688          -0.74120
## 1500 meter buffer          0.012890       -0.8887          -0.10820
## 1750 meter buffer          0.024470       -0.9706           0.07963
## 2500 meter buffer          0.207500       -0.8769           0.32580
## 2250 meter buffer          0.257900       -0.8396           0.10360
## 4250 meter buffer         -0.498100       -0.8505           0.68850
## 2000 meter buffer          0.080410       -0.9469           0.09995
## 750 meter buffer          -0.027740       -0.6650          -0.73090
## 4000 meter buffer         -0.264100       -0.5710           0.59510
## 2750 meter buffer          0.029550       -0.7075           0.32390
## 3750 meter buffer         -0.155200       -0.3596           0.42950
## 500 meter buffer          -0.182500       -0.6500          -0.39050
## 3000 meter buffer          0.069910       -0.5419           0.33260
## 250 meter buffer          -0.056660       -0.9442          -0.44440
## 3500 meter buffer         -0.089060       -0.3220           0.28860
## 3250 meter buffer          0.111800       -0.3295           0.29970
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik  AICc delta
## 4750 meter buffer           2.11400       -2.1910 18 -149.588 338.4  0.00
## 5000 meter buffer           2.24100       -1.7140 18 -151.928 343.1  4.68
## 4500 meter buffer           1.65400       -1.8260 18 -152.395 344.0  5.62
## 1000 meter buffer           0.32020       -0.5817 18 -155.295 349.8 11.41
## 1250 meter buffer           0.47220       -0.7001 18 -155.337 349.9 11.50
## 1500 meter buffer           0.72530       -1.0730 18 -155.370 350.0 11.56
## 1750 meter buffer           0.68240       -1.1560 18 -155.922 351.1 12.67
## 2500 meter buffer           0.71000       -1.3560 18 -156.600 352.4 14.02
## 2250 meter buffer           0.62500       -1.2220 18 -157.202 353.6 15.23
## 4250 meter buffer           1.11700       -0.8839 18 -158.164 355.5 17.15
## 2000 meter buffer           0.65430       -1.0220 18 -158.480 356.2 17.78
## 750 meter buffer            0.38590       -0.5004 18 -159.420 358.1 19.67
## 4000 meter buffer           0.77010       -0.5758 18 -161.281 361.8 23.39
## 2750 meter buffer           0.75080       -1.0640 18 -161.352 361.9 23.53
## 3750 meter buffer           0.85350       -0.6662 18 -162.517 364.2 25.86
## 500 meter buffer            0.37150        0.1016 18 -162.994 365.2 26.81
## 3000 meter buffer           0.69030       -0.9685 18 -163.140 365.5 27.10
## 250 meter buffer           -0.00829        0.1530 18 -163.822 366.9 28.47
## 3500 meter buffer           0.85220       -0.7915 18 -165.724 370.7 32.27
## 3250 meter buffer           0.75900       -0.4774 18 -166.878 373.0 34.58
##                   weight
## 4750 meter buffer  0.855
## 5000 meter buffer  0.082
## 4500 meter buffer  0.052
## 1000 meter buffer  0.003
## 1250 meter buffer  0.003
## 1500 meter buffer  0.003
## 1750 meter buffer  0.002
## 2500 meter buffer  0.001
## 2250 meter buffer  0.000
## 4250 meter buffer  0.000
## 2000 meter buffer  0.000
## 750 meter buffer   0.000
## 4000 meter buffer  0.000
## 2750 meter buffer  0.000
## 3750 meter buffer  0.000
## 500 meter buffer   0.000
## 3000 meter buffer  0.000
## 250 meter buffer   0.000
## 3500 meter buffer  0.000
## 3250 meter buffer  0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For red fox the top model is 4750 m buffer

Model summary 4750m

Let’s get the model summary

summary(fox_mods$`4750 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(red_fox, absent_red_fox) ~ scale(harvest) + scale(pipeline) +  
##     scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##    335.2    397.2   -149.6    299.2      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 10.88    3.298   
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -5.48167    1.42811  -3.838 0.000124 ***
## scale(harvest)             0.33898    0.36242   0.935 0.349620    
## scale(pipeline)            0.55450    0.58497   0.948 0.343173    
## scale(roads)               1.77163    0.62648   2.828 0.004685 ** 
## scale(seismic_lines)      -0.77875    0.35358  -2.203 0.027629 *  
## scale(seismic_lines_3D)   -0.25981    0.62085  -0.418 0.675594    
## scale(trails)             -1.12241    0.52768  -2.127 0.033413 *  
## scale(transmission_lines)  0.51430    0.28582   1.799 0.071957 .  
## scale(veg_edges)           2.11425    0.61075   3.462 0.000537 ***
## scale(wells)              -2.19119    0.74700  -2.933 0.003354 ** 
## scale(osm_industrial)      0.06849    0.50631   0.135 0.892402    
## scale(lc_grassland)        0.45871    0.31291   1.466 0.142671    
## scale(lc_coniferous)      -1.73437    0.97937  -1.771 0.076576 .  
## scale(lc_broadleaf)       -0.59669    0.71105  -0.839 0.401381    
## scale(lc_mixed)           -2.13469    0.56671  -3.767 0.000165 ***
## scale(lc_developed)       -1.94622    0.69304  -2.808 0.004981 ** 
## scale(lc_shrub)           -0.21928    0.60779  -0.361 0.718262    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looks good so far

Subset Models

Add later if deemed necessary.

Top model

White-tailed deer

Global models

deer_mods <- osm_final_df_2021_2022 %>%
  
  # use purrr map to make global models for all other buffer sizes
  purrr::map(
    ~.x %>%

      # have to include the `` around the white-tailed_deer or R won't recognize it as a variable because of the -
     glmmTMB::glmmTMB(cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~
                        
                        # HFI
                        scale(harvest) +
                        scale(pipeline) +
                        scale(roads) +
                        scale(seismic_lines) +
                        scale(seismic_lines_3D) +
                        scale(trails) +
                        scale(transmission_lines) +
                        scale(veg_edges) +
                        scale(wells) +
                        scale(osm_industrial) +
                        
                        # VEG covariates in numerical order
                        scale(lc_grassland) +
                        scale(lc_coniferous) +
                        scale(lc_broadleaf) +
                        scale(lc_mixed) +
                        scale(lc_developed) +
                        scale(lc_shrub) +
                        
                        # Random effect of array
                        (1|array),
                   data = .,
                   family = 'binomial'))


# model selection
model.sel(deer_mods)
## Model selection table 
##                   cnd((Int)) dsp((Int)) cnd(scl(hrv)) cnd(scl(lc_brd))
## 1500 meter buffer    -0.2043          +      -0.14420         0.145500
## 1250 meter buffer    -0.1948          +      -0.15150         0.173400
## 1750 meter buffer    -0.2032          +      -0.15680         0.113000
## 4000 meter buffer    -0.2110          +      -0.29720        -0.072220
## 3750 meter buffer    -0.2075          +      -0.28380        -0.049260
## 3500 meter buffer    -0.2083          +      -0.26840        -0.018200
## 4250 meter buffer    -0.2113          +      -0.30650        -0.086260
## 2000 meter buffer    -0.2051          +      -0.17040         0.139300
## 4500 meter buffer    -0.2117          +      -0.30190        -0.139500
## 3250 meter buffer    -0.2060          +      -0.23600        -0.028290
## 4750 meter buffer    -0.2111          +      -0.29520        -0.199600
## 5000 meter buffer    -0.2104          +      -0.29120        -0.255600
## 2250 meter buffer    -0.2013          +      -0.19010         0.116300
## 3000 meter buffer    -0.2015          +      -0.20990        -0.001321
## 2750 meter buffer    -0.2006          +      -0.19840         0.061800
## 2500 meter buffer    -0.2004          +      -0.18860         0.125000
## 1000 meter buffer    -0.1789          +      -0.16970         0.237800
## 750 meter buffer     -0.1827          +      -0.12340         0.159000
## 500 meter buffer     -0.1890          +      -0.04165        -0.040670
## 250 meter buffer     -0.1830          +       0.02303        -0.103900
##                   cnd(scl(lc_cnf)) cnd(scl(lc_dvl)) cnd(scl(lc_grs))
## 1500 meter buffer         -0.17350        -0.007861          0.15040
## 1250 meter buffer         -0.12080         0.038340          0.13320
## 1750 meter buffer         -0.20820        -0.026170          0.14180
## 4000 meter buffer         -0.52070        -0.136700          0.28260
## 3750 meter buffer         -0.48300        -0.157200          0.24990
## 3500 meter buffer         -0.41060        -0.117800          0.24370
## 4250 meter buffer         -0.56400        -0.179700          0.29920
## 2000 meter buffer         -0.14010         0.001791          0.17340
## 4500 meter buffer         -0.65010        -0.212100          0.31300
## 3250 meter buffer         -0.38630        -0.079840          0.23610
## 4750 meter buffer         -0.73840        -0.208000          0.30060
## 5000 meter buffer         -0.83550        -0.215400          0.27990
## 2250 meter buffer         -0.16340        -0.025840          0.17270
## 3000 meter buffer         -0.33820        -0.055770          0.22780
## 2750 meter buffer         -0.23950        -0.027520          0.20890
## 2500 meter buffer         -0.16090        -0.037020          0.20250
## 1000 meter buffer         -0.02351         0.100100          0.10230
## 750 meter buffer          -0.08746         0.155200          0.06230
## 500 meter buffer          -0.37420         0.078840         -0.01527
## 250 meter buffer          -0.56080        -0.001368         -0.08770
##                   cnd(scl(lc_mxd)) cnd(scl(lc_shr)) cnd(scl(osm_ind))
## 1500 meter buffer          0.28160         0.033630            0.4141
## 1250 meter buffer          0.28580         0.054670            0.4178
## 1750 meter buffer          0.26060         0.006141            0.3574
## 4000 meter buffer          0.21600        -0.118900            0.1799
## 3750 meter buffer          0.23940        -0.095690            0.1647
## 3500 meter buffer          0.26370        -0.061460            0.1660
## 4250 meter buffer          0.16210        -0.139000            0.1654
## 2000 meter buffer          0.30840         0.053940            0.3139
## 4500 meter buffer          0.09763        -0.193100            0.1473
## 3250 meter buffer          0.23790        -0.056530            0.1791
## 4750 meter buffer          0.02768        -0.264300            0.1247
## 5000 meter buffer         -0.04817        -0.329700            0.1121
## 2250 meter buffer          0.28140         0.037560            0.2715
## 3000 meter buffer          0.24280        -0.036920            0.2004
## 2750 meter buffer          0.26830         0.010970            0.2197
## 2500 meter buffer          0.27030         0.047620            0.2355
## 1000 meter buffer          0.29890         0.066870            0.3539
## 750 meter buffer           0.16780         0.014890            0.3334
## 500 meter buffer          -0.04650        -0.137800            0.2630
## 250 meter buffer          -0.16920        -0.236400            0.1880
##                   cnd(scl(ppl)) cnd(scl(rds)) cnd(scl(ssm_lns_3D))
## 1500 meter buffer      0.058750    -0.1019000              -0.4396
## 1250 meter buffer      0.095410    -0.0485900              -0.4204
## 1750 meter buffer      0.014140    -0.0007351              -0.3973
## 4000 meter buffer      0.002368    -0.1157000              -0.2863
## 3750 meter buffer      0.101400    -0.0711400              -0.3089
## 3500 meter buffer      0.132100    -0.0364700              -0.3089
## 4250 meter buffer     -0.020040    -0.0965900              -0.2656
## 2000 meter buffer     -0.005024     0.0196300              -0.3705
## 4500 meter buffer     -0.089730    -0.0727300              -0.2495
## 3250 meter buffer      0.126100    -0.0133700              -0.2967
## 4750 meter buffer     -0.195000    -0.0433200              -0.2262
## 5000 meter buffer     -0.279500    -0.0379400              -0.2021
## 2250 meter buffer     -0.007749     0.0586200              -0.3343
## 3000 meter buffer      0.057480     0.0145900              -0.2858
## 2750 meter buffer      0.020450     0.0491900              -0.2936
## 2500 meter buffer      0.018310     0.0742400              -0.3052
## 1000 meter buffer      0.129300     0.0047080              -0.3375
## 750 meter buffer       0.114200    -0.0349300              -0.2754
## 500 meter buffer       0.096410    -0.0344700              -0.2306
## 250 meter buffer       0.105900    -0.0198700              -0.2150
##                   cnd(scl(ssm_lns)) cnd(scl(trl)) cnd(scl(trn_lns))
## 1500 meter buffer          -0.14640       0.11130         -0.053540
## 1250 meter buffer          -0.13820       0.10560         -0.104600
## 1750 meter buffer          -0.14340       0.11520         -0.003427
## 4000 meter buffer          -0.08145       0.20080         -0.126600
## 3750 meter buffer          -0.07416       0.20910         -0.181200
## 3500 meter buffer          -0.06994       0.19850         -0.225800
## 4250 meter buffer          -0.10260       0.19180         -0.135000
## 2000 meter buffer          -0.13400       0.11600         -0.012350
## 4500 meter buffer          -0.09499       0.18090         -0.102900
## 3250 meter buffer          -0.04894       0.18850         -0.205200
## 4750 meter buffer          -0.07050       0.16710         -0.061650
## 5000 meter buffer          -0.04943       0.15980         -0.044950
## 2250 meter buffer          -0.09048       0.13470         -0.041230
## 3000 meter buffer          -0.03858       0.16090         -0.128600
## 2750 meter buffer          -0.05687       0.15650         -0.081490
## 2500 meter buffer          -0.07414       0.15470         -0.080920
## 1000 meter buffer          -0.12240       0.12710         -0.124100
## 750 meter buffer           -0.04965       0.06837          0.008704
## 500 meter buffer            0.01582       0.08367          0.071480
## 250 meter buffer           -0.05631       0.16840         -0.011320
##                   cnd(scl(veg_edg)) cnd(scl(wll)) df   logLik   AICc delta
## 1500 meter buffer         -0.242100       0.39620 18 -537.053 1113.3  0.00
## 1250 meter buffer         -0.249000       0.32870 18 -541.160 1121.5  8.22
## 1750 meter buffer         -0.278300       0.37420 18 -541.829 1122.9  9.55
## 4000 meter buffer          0.001398       0.40600 18 -544.404 1128.0 14.70
## 3750 meter buffer          0.016500       0.36800 18 -544.421 1128.1 14.74
## 3500 meter buffer         -0.021160       0.34410 18 -545.857 1130.9 17.61
## 4250 meter buffer          0.060980       0.38260 18 -545.869 1130.9 17.63
## 2000 meter buffer         -0.274900       0.37520 18 -546.342 1131.9 18.58
## 4500 meter buffer          0.104900       0.38530 18 -546.625 1132.5 19.14
## 3250 meter buffer         -0.085550       0.35550 18 -548.339 1135.9 22.57
## 4750 meter buffer          0.124800       0.39700 18 -548.938 1137.1 23.77
## 5000 meter buffer          0.138100       0.43710 18 -549.440 1138.1 24.77
## 2250 meter buffer         -0.232100       0.39110 18 -550.150 1139.5 26.19
## 3000 meter buffer         -0.133500       0.36960 18 -551.018 1141.2 27.93
## 2750 meter buffer         -0.194200       0.38230 18 -551.979 1143.2 29.85
## 2500 meter buffer         -0.183800       0.35610 18 -552.560 1144.3 31.01
## 1000 meter buffer         -0.199000       0.27410 18 -552.713 1144.6 31.32
## 750 meter buffer          -0.239500       0.25260 18 -568.054 1175.3 62.00
## 500 meter buffer          -0.221400       0.21880 18 -573.059 1185.3 72.01
## 250 meter buffer          -0.215200       0.05339 18 -574.675 1188.6 75.24
##                   weight
## 1500 meter buffer  0.974
## 1250 meter buffer  0.016
## 1750 meter buffer  0.008
## 4000 meter buffer  0.001
## 3750 meter buffer  0.001
## 3500 meter buffer  0.000
## 4250 meter buffer  0.000
## 2000 meter buffer  0.000
## 4500 meter buffer  0.000
## 3250 meter buffer  0.000
## 4750 meter buffer  0.000
## 5000 meter buffer  0.000
## 2250 meter buffer  0.000
## 3000 meter buffer  0.000
## 2750 meter buffer  0.000
## 2500 meter buffer  0.000
## 1000 meter buffer  0.000
## 750 meter buffer   0.000
## 500 meter buffer   0.000
## 250 meter buffer   0.000
## Models ranked by AICc(x) 
## Random terms (all models): 
##   cond(1 | array)

For deer the top model was also the 1500 buffer

Model summary 1500m

Let’s get the model summary

summary(deer_mods$`1500 meter buffer`)
##  Family: binomial  ( logit )
## Formula:          
## cbind(`white-tailed_deer`, `absent_white-tailed_deer`) ~ scale(harvest) +  
##     scale(pipeline) + scale(roads) + scale(seismic_lines) + scale(seismic_lines_3D) +  
##     scale(trails) + scale(transmission_lines) + scale(veg_edges) +  
##     scale(wells) + scale(osm_industrial) + scale(lc_grassland) +  
##     scale(lc_coniferous) + scale(lc_broadleaf) + scale(lc_mixed) +  
##     scale(lc_developed) + scale(lc_shrub) + (1 | array)
## Data: .
## 
##      AIC      BIC   logLik deviance df.resid 
##   1110.1   1172.1   -537.1   1074.1      214 
## 
## Random effects:
## 
## Conditional model:
##  Groups Name        Variance Std.Dev.
##  array  (Intercept) 1.107    1.052   
## Number of obs: 232, groups:  array, 6
## 
## Conditional model:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               -0.204311   0.432747  -0.472   0.6368    
## scale(harvest)            -0.144219   0.068902  -2.093   0.0363 *  
## scale(pipeline)            0.058746   0.107685   0.546   0.5854    
## scale(roads)              -0.101872   0.146611  -0.695   0.4872    
## scale(seismic_lines)      -0.146367   0.069900  -2.094   0.0363 *  
## scale(seismic_lines_3D)   -0.439551   0.096937  -4.534 5.78e-06 ***
## scale(trails)              0.111275   0.053479   2.081   0.0375 *  
## scale(transmission_lines) -0.053540   0.118825  -0.451   0.6523    
## scale(veg_edges)          -0.242143   0.130037  -1.862   0.0626 .  
## scale(wells)               0.396152   0.088212   4.491 7.09e-06 ***
## scale(osm_industrial)      0.414131   0.068684   6.030 1.64e-09 ***
## scale(lc_grassland)        0.150428   0.076396   1.969   0.0489 *  
## scale(lc_coniferous)      -0.173514   0.238731  -0.727   0.4673    
## scale(lc_broadleaf)        0.145490   0.169488   0.858   0.3907    
## scale(lc_mixed)            0.281594   0.155049   1.816   0.0693 .  
## scale(lc_developed)       -0.007861   0.118023  -0.067   0.9469    
## scale(lc_shrub)            0.033628   0.139233   0.242   0.8091    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Looks okay for now

Subset models

Add later if deemed necessary

Top model

Top buffers

For convenience here is a list of the top buffer widths for each species

  • Black bear - 250m
  • Caribou - 3000m
  • Coyote - 2250m
  • Fisher - 1000m
  • Grey wolf - 500m
  • Lynx - 1250m
  • Moose - 250m
  • Red fox - 4750
  • White-tailed deer - 1500m

Plots for all species

Buffer size

What I want to do here is create a figure that plots a point for each species at the buffer size that was in the top model.

Data frame

Before we can do this I need to make a data frame that include the species, top buffer width, and any other info I may want in the plot.

buffer_data <- tibble(
  
  #define number of rows
  .rows = 9,
  
  # create species column
  species = c('Black bear',
              'Caribou',
              'Coyote',
              'Fisher',
              'Grey wolf',
              'Lynx',
              'Moose',
              'Red fox',
              'White-tailed deer'),
  
  # create buffer column w/ buffer corresponding to the top model for each species
  buffer_width = c(250,
                   3000,
                   2250,
                   1000,
                   500,
                   1250,
                   250,
                   4750,
                   1500),
  
   # add phylopic uuid for each species for plotting 
  image = c(get_uuid(name = 'Ursus americanus'),
           get_uuid(name = 'Rangifer tarandus'),
           get_uuid(name = 'Canis latrans'),
           get_uuid(name = 'Pekania pennanti'),
           get_uuid(name = 'Canis lupus'),
           get_uuid(name = 'Lynx lynx'),
           get_uuid(name = 'Alces alces'),
           get_uuid(name = 'Vulpes vulpes'),
           get_uuid(name ='Odocoileus virginianus'))
  
)

Plot

ggplot(buffer_data,
       aes( x = species,
            y = buffer_width)) +
  
  # add points as phylopic silhouettes
  geom_phylopic(aes(uuid = image),
                size = 500) +
  
  scale_y_continuous(breaks = seq(0, 5000, by = 250)) +
  
  coord_cartesian(ylim = c(0, 4900))

This looks alright but I don’t love all the default phylopic images that are used.

Fix phylopics

We can go to the Phylopic website and check out other options and then use the code below to search for the uuid from those options for the ones we want to replace

# Caribou could be bigger
get_uuid(name = 'Rangifer tarandus',
         n = 6) 
## [1] "09c5d6bf-8ed9-4dc3-86f1-6b75dbe63952"
## [2] "a7e8f29e-a4e2-499a-b913-491f78f44f1d"
## [3] "ad1de77a-c89b-41ad-801c-1ef47c960c83"
## [4] "d2e268c1-1937-4308-9fa4-9adc55cc5e8b"
## [5] "e6e864fd-8e3d-435f-9db3-dc6869c589f1"
## [6] "19864ae2-b740-4a41-9f91-e3a195ee7cbb"
# I like the 5th one best

# Coyote doesn't need to have a food item
get_uuid(name = 'Canis latrans',
         n = 6) 
## [1] "113d2520-9f92-456f-b305-52ee3986172d"
## [2] "3492f4ca-01f0-4609-a84f-084a84bf4e95"
## [3] "5a0398e3-a455-4ca6-ba86-cf3f1b25977a"
## [4] "76c8fdec-d0af-47b9-b949-dc610419d832"
## [5] "da5faa63-085f-4523-a542-e71cb386c999"
## [6] "e6a2fa4b-85df-43b4-989c-34a65ba7eee3"
# I like the 6th one best

# Grey wolf shouldn't be the same as coyote
get_uuid(name = 'Canis lupus',
         n = 5) 
## [1] "113d2520-9f92-456f-b305-52ee3986172d"
## [2] "11658f8c-e0c2-4612-85ef-bdc44acdae0b"
## [3] "3492f4ca-01f0-4609-a84f-084a84bf4e95"
## [4] "5a0398e3-a455-4ca6-ba86-cf3f1b25977a"
## [5] "76c8fdec-d0af-47b9-b949-dc610419d832"
# I like 2 or 5

# Moose needs a body
get_uuid(name = 'Alces alces',
         n = 3) 
## [1] "df2d0ad0-adb0-49d7-afe5-edc6cad21064"
## [2] "1a20a65d-1342-4833-a9dd-1611b9fb383c"
## [3] "74eab34a-498c-4614-aece-f02361874f79"
# 3 is best

# red fox is too big
get_uuid(name = 'Vulpes vulpes',
         n = 7) 
## [1] "76352962-1eeb-4197-acdd-e3c7eeab839d"
## [2] "9b3e3567-2238-40bf-ab5d-bcb6ce7b32d5"
## [3] "a5e2a085-c895-4fdd-b39e-bbac8ca94d7d"
## [4] "d67d3bf6-3509-4ab6-819a-cd409985347e"
## [5] "40541f06-01d8-4585-a271-40268d24057d"
## [6] "40a5954e-ede1-4af6-9f65-209d1780e602"
## [7] "51b1b6e4-129d-41a6-bbbd-c3fab459c25f"
# 4 is best

# deer is weird
get_uuid(name = 'Odocoileus virginianus',
         n = 6) 
## [1] "4584be20-4514-4673-a3e8-97e2a6a10e57"
## [2] "49a5a5db-047e-4e17-849b-9f96a93f0d2b"
## [3] "56f6fdb2-15d0-43b5-b13f-714f2cb0f5d0"
## [4] "9df2ed00-b918-4d3a-82bd-1a94283f7c0f"
## [5] "6038e80c-398d-47b2-9a69-2b9edf436f64"
## [6] "8569838c-c725-4772-b0a3-b5eb04baaada"
# 5 is best

Let’s alter the code from above to change the phylopic images

buffer_data <- tibble(
  
  #define number of rows
  .rows = 9,
  
  # create species column
  species = c('Black bear',
              'Caribou',
              'Coyote',
              'Fisher',
              'Grey wolf',
              'Lynx',
              'Moose',
              'Red fox',
              'White-tailed deer'),
  
  # create buffer column w/ buffer corresponding to the top model for each species
  buffer_width = c(250,
                   3000,
                   2250,
                   1000,
                   500,
                   1250,
                   250,
                   4750,
                   1500),
  
   # add phylopic uuid for each species for plotting 
  image = c(get_uuid(name = 'Ursus americanus'),
            'e6e864fd-8e3d-435f-9db3-dc6869c589f1', # caribou 5
            'e6a2fa4b-85df-43b4-989c-34a65ba7eee3', # coyote 6
            get_uuid(name = 'Pekania pennanti'),
            '76c8fdec-d0af-47b9-b949-dc610419d832', # wolf 5
            get_uuid(name = 'Lynx lynx'),
            '74eab34a-498c-4614-aece-f02361874f79', # moose 3
            'd67d3bf6-3509-4ab6-819a-cd409985347e', # red fox 4
            '6038e80c-398d-47b2-9a69-2b9edf436f64') # white-tailed deer 6
  
)

Plot again

Using the code from above with the fixed images let’s plot it again

buffer_plot <- ggplot(buffer_data,
                      aes( x = species,
                           y = buffer_width)) +
  
  # add points as phylopic silhouettes
  geom_phylopic(aes(uuid = image),
                size = 400) +
  
  # set axis limits and breaks
  scale_y_continuous(breaks = seq(0, 5000, by = 500)) +
  
  # wrap long species names
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  
  # set graph limits
  coord_cartesian(ylim = c(0,4900)) +
  
  # make nicer labels
  labs(y = 'Buffer width (meters)') +
  
  # set theme
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

buffer_plot 

BEAUTIFUL!!! If I do say so myself

Save plot

ggsave('figures/OSM_species_buffers_2022.jpeg',
       buffer_plot,
       width = 18,
       height = 12,
       units = 'in',
       dpi = 600)

Body size plot

Now I want to see if their are trends between body size and which bfuffer distance performed best. I’ll first need to add a body size column to the data.

Mutate data

To do this I’m going to google the avergae body size of each species and then add it as a column using mutate

** IMPORTANT - this is just for data visualization purposes!. Since we don’t actually have data on body size of the individuals in our study I can’t run any analysis on this but I can use it to visualize if there appears to be any trends worth investigating.

# first add body size to buffer data
body_size_data <- buffer_data %>% 
  
  # use mutate to create new column
  mutate(body_size = case_when(species == 'Black bear' ~ 181,
                               species == 'Caribou' ~ 355,
                               species == 'Coyote' ~ 15,
                               species == 'Fisher' ~ 3.5,
                               species == 'Grey wolf' ~ 39,
                               species == 'Lynx' ~ 11,
                               species == 'Moose' ~ 555,
                               species == 'Red fox' ~ 5.2,
                               species == 'White-tailed deer' ~ 200),
         log_body_size = log(body_size))

Plot

Now let’s plot this as a scatterplot

body_size_plot <- ggplot(body_size_data, aes(x = log_body_size,
                                             y = buffer_width)) +
  
  # add points as phylopic silhouettes
  geom_phylopic(aes(uuid = image),
                size = 400) +
  
  labs(x = 'Log body size (kg)',
       y = 'Buffer width (meters)') +
  
  # set theme
  theme_bw() +
  theme(panel.grid.minor = element_blank(),
        axis.text = element_text(size = 14),
        axis.title = element_text(size = 16))

body_size_plot

Save plot

ggsave('figures/OSM_body_size_buffers_2022.jpg',
       body_size_plot,
       width = 18,
       height = 12,
       units = 'in',
       dpi = 600)

Odds ratio plots

What we need to do now is extract the coefficient estimates from each top model, as well as the 95% CI so we can plot them for easier visualization and understanding of the data

Extract odds ratios

The confint() function will extract the coefficients and CI intervals from a model, so what we need to do is make a list of all the models, then use the map() function in purrr to apply the confint() function to all the models and get the coefficients. We want this to result in a tibble that has a column for the HFI feature (we aren’t plotting the lc_class data for this report), the upper and lower CI, and the coefficient estimate.

In order to do this we have to do a bit of data wrangling, currently this isn’t the most pleasing way to accomplish the desired outcome, but it works.

# This is also a dog shit way to do this but I need to get this done

# make a list of the top models for each species
top_models <- list(black_bear_mods$`250 meter buffer`,
                   caribou_mods$`3000 meter buffer`,
                   coyote_mods$`2250 meter buffer`,
                   fisher_mods$`1000 meter buffer`,
                   wolf_mods$`500 meter buffer`,
                   lynx_mods$`1250 meter buffer`,
                   moose_mods$`250 meter buffer`,
                   fox_mods$`4750 meter buffer`,
                   deer_mods$`1500 meter buffer`) %>%  
  
  # pipe into purrr to create coefficient table for all models
  purrr::map(
    ~.x %>% 
      
      # extract the coefficients and upper and lower CI
      confint() %>% 
      
      # format resulting object as a tibble data frame
      as_tibble() %>%
      
      # subset to just the HFI variables for these plots
      slice_head(n = 11) %>% 
      
      # remove first row which is the intercept
      slice_tail(n = 10) %>% 
      
      # add a column where we can put the feature names
      rowid_to_column() %>% 
      
      # rename the columns for plotting
      rename('lower' = `2.5 %`,
             'upper' = `97.5 %`,
             'estimate' = Estimate,
             'feature' = rowid) %>% 
      
      # rename the entries to features, need to look at the order the features are in from the model summary and ensure it matches
      mutate(feature = as.factor(feature),
             feature = recode(feature,
                              '1' = 'harvest',
                              '2' = 'pipeline',
                              '3' = 'roads',
                              '4' = 'seismic_lines',
                              '5' = 'seismic_lines_3D',
                              '6' = 'trails',
                              '7' = 'transmission_lines',
                              '8' = 'veg_edges',
                              '9' = 'wells',
                              '10' = 'osm_industrial'))) %>% 
  
  # set the names of each resulting tibble data frame to the species name
  purrr::set_names('Black bear',
                   'Caribou',
                   'Coyote',
                   'Fisher',
                   'Grey wolf',
                   'Lynx',
                   'Moose',
                   'Red fox',
                   'White-tailed deer')

Merge data frames

Now we have a data frame with the odds ratios for each species, but if we want these on a plot together we need them all in one data frame.

To merge data into one data frame we can use the list_rbind() function from the purrr package which will take each element of the list and stack them on top of one another just like rbind does with data frames, and if we use the names_to argument we can extract the names of the list elements and assign them to a column so we know which data comes from which species model (list element)

In this code I also add a new column called uuid which contains the image id (uuid) for a phylopic silhouette of each species that I may want to use for plotting

Phylopic.or is an open source online database of silhouettes various contributors have created for use. There is an R package that works with this data called rphylopic; we can use the get_uuid() function from this package to extract the data for a silhouette for each species we want, which is what I’ve done here.

# combine all list elements 
coeffs_df_all <- list_rbind(top_models,
                            names_to = 'species') %>% 
  
  # change species to a factor for plotting
  mutate(species = as.factor(species),
         
         # add phylopic uuid for each species for plotting 
         # the uuid is extracted using getuuid with the species name as name = ''
         uuid = case_when(species == 'Black bear' ~ get_uuid(name = 'Ursus americanus'),
                          species == 'Caribou' ~ get_uuid(name = 'Rangifer tarandus'),
                          species == 'Coyote' ~ get_uuid(name = 'Canis latrans'),
                          species == 'Fisher' ~ get_uuid(name = 'Pekania pennanti'),
                          species == 'Grey wolf' ~ get_uuid(name = 'Canis lupus'),
                          species == 'Lynx' ~ get_uuid(name = 'Lynx lynx'),
                          species == 'Moose' ~ get_uuid(name = 'Alces alces'),
                          species == 'Red fox' ~ get_uuid(name = 'Vulpes vulpes'),
                          species == 'White-taield deer' ~ get_uuid(name ='Odocoileus virginianus'))) 

Plot odds all species

Now let’s explore some different options to plot the coefficients

geom_phylopic

Let’s try plotting all the species on one plot using ggplot()

# provide data and mapping aesthetics
ggplot(coeffs_df_all, aes(x = feature, 
                          y = estimate, 
                          group = uuid)) +
  
   geom_errorbar(aes(ymin = lower, 
                     ymax = upper, 
                     color = feature),
                width = 0.4,
                linewidth = 0.5,
                position = position_dodge(width = 0.9)) +
  
    # add points for each estimate for each covariate and use position = position_dodge to shift the points so all the species don't plot on top of one another
  geom_phylopic(aes(x = feature, 
                 y = estimate,
                 uuid = uuid),
             position = position_dodge(width = 0.9),
             size = 0.3)
## Warning: Removed 10 rows containing missing values.

This works but unfortunately is very messy with the images as the points since they don’t all have the same centerline and are different sizes.

I try with the code below to select better silhouettes

# combine all list elements 
coeffs_df_all <- list_rbind(top_models,
                            names_to = 'species') %>% 
  
  # change species to a factor for plotting
  mutate(species = as.factor(species),
         
         # add phylopic uuid for each species for plotting 
         # the uuid is extracted using getuuid with the species name as name = ''
         uuid = case_when(species == 'Black bear' ~ get_uuid(name = 'Ursus americanus'),
                           species == 'Caribou' ~ 'e6e864fd-8e3d-435f-9db3-dc6869c589f1', # caribou 5
                          species == 'Coyote' ~ 'e6a2fa4b-85df-43b4-989c-34a65ba7eee3', # coyote 6
                          species == 'Fisher' ~ '735066c6-2f3e-4f97-acb1-06f55ae075c9',
                          species == 'Grey wolf' ~ '76c8fdec-d0af-47b9-b949-dc610419d832', # wolf 5
                          species == 'Lynx' ~ '74eab34a-498c-4614-aece-f02361874f79', # moose 3
                          species == 'Moose' ~ '74eab34a-498c-4614-aece-f02361874f79',
                          species == 'Red fox' ~ 'd67d3bf6-3509-4ab6-819a-cd409985347e', # red fox 4
                          species == 'White-taield deer' ~ '6038e80c-398d-47b2-9a69-2b9edf436f64')) # white-tailed deer 6

Try plotting all again

ggplot(coeffs_df_all, aes(x = feature, 
                          y = estimate, 
                          group = uuid)) +
  
   geom_errorbar(aes(ymin = lower, 
                     ymax = upper, 
                     color = feature),
                width = 0.4,
                linewidth = 0.5,
                position = position_dodge(width = 1.2)) +
  
    # add points for each estimate for each covariate and use position = position_dodge to shift the points so all the species don't plot on top of one another
  geom_phylopic(aes(x = feature, 
                 y = estimate,
                 uuid = uuid),
             position = position_dodge(width = 1.2),
             size = 0.2)
## Warning: Removed 10 rows containing missing values.
## Warning: `position_dodge()` requires non-overlapping x intervals

Unfortunately they still don’t center on the lines like they should so will try a different plot

Facet_wrap

ggplot(coeffs_df_all,
       aes(x = feature,
           y = estimate)) +
  
  # add points
  geom_point() +
  
  # add errorbar
  geom_errorbar(aes(ymin = lower, 
                     ymax = upper),
                width = 0.4,
                linewidth = 0.5) +
  
  facet_wrap(vars(species)) +
  
  # wrap long axis names
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  
  # adjust theme elements
  theme(axis.text.x = element_text(angle = 90,
                                   hjust = 1),
        axis.ticks.x = element_blank())

Plot_model

If all else fails can use plot_model function

Example with black bear model
library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
plot_model(black_bear_mods$`250 meter buffer`,
           vline.color = 'black')